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ABSTRACT 

Aims. We present a panchromatic study, involving a multiple technique approach, of the circumstellar disc surrounding 
the T Tauri star IM Lupi (Sz 82). 

Methods. We have undertaken a comprehensive observational study of IM Lupi using photometry, spectroscopy, millime- 
tre interferometry and multi-wavelength imaging. For the first time, the disc is resolved from optical and near-infrared 
wavelengths in scattered light, to the millimetre regime in thermal emission. Our data-set, in conjunction with existing 
photometric data, provides an extensive coverage of the spectral energy distribution, including a detailed spectrum of 
the silicate emission bands. We have performed a simultaneous modelling of the various observations, using the radia- 
tive transfer code MCFOST, and analysed a grid of models over a large fraction of the parameter space via Bayesian 
inference. 

Results. We have constructed a model that can reproduce all of the observations of the disc. Our analysis illustrates the 
importance of combining a wide range of observations in order to fully constrain the disc model, with each observation 
providing a strong constraint only on some aspects of the disc structure and dust content. Quantitative evidence of 
dust evolution in the disc is obtained: grain growth up to millimetre-sized particles, vertical stratification of dust grains 
with micrometric grains close to the disc surface and larger grains which have settled towards the disc midplane, and 
possibly the formation of fluffy aggregates and/or ice mantles around grains. 

Key words, circumstellar matter — Accretion discs — planetary systems: proto-planetary discs — Radiative transfer 
— stars: formation, individual: IM Lupi 



1. Introduction 

During the first stages of planet fo rmation following the 
core nucleated accretion scenario ()Lissauer fc Stevenson! 
|2007| ). evolution of dust grains within the protoplanetary 
disc surrounding the central forming object is expected. 
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Micrometre size dust grains will start to grow by co- 
agulation during low relative velocity collisions, leading 
to the formation of l a rger, potentially fluf f y aggregates 
(iBeckwith et all 120001: iDominik et~aH l2007t iNatta et all 
120071 and references therein) that will ultimately give birth 
to kilometre size planetesimals promoting gravitational fo- 
cusing and eventually leading to planets. Timescales for 
the formation of these aggregates strongly depend on the 
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Fig. 2. NICMOS F160W image in log stretch. The field of 
view is If xfl" with North up and East to the left. The 
dashed box corresponds to the field of view in Fig. [T] The 
central white circle represents the 0.3" radius coronagraphic 
obscuration. The full green lines indicate the upper and 
lower scattering surfaces of the disc (separated by the dark 
lane which corresponds to the disc midplane). The green 
dashed circle represents the possible envelope surrounding 
the disc. 



physics of aggregation as well as the physical conditions 
within the disc, such as differential velocities between grains 
and grain-gas interactions. 

In parallel to grain growth, dust grains are expected 
to settle towards the disc midplane and then to migrate 
inward as a result of the c onjugate actions of the stel- 
lar gravity and gas drag (e.g . iBarriere-Fouchet et al.|[2005t 
iFromang fc Papaloizoul2006r ). This settling is highly depen- 
dent on the grain size. Small grains (< I /zm) are strongly 
coupled to the gas, they follow its motion, and do not set- 
tle at all. Conversely, grains in the mm-cm regime are con- 
siderably slowed down by the gas drag. They completely 
decouple and settle very rapidly in a thin midplane. 

The result is the formatio n of a dust sub-disc close 
to the equator i al pla ne (e.g.. ISafronov fc Zviaginal Il969t 
iDubrulle et all II995D that may become unstable and 
produce planetesimals when t he density of the dust 
layer exceeds the gas one (e.g.. iGoldreich fc Ward! fl973t 
ISchrapler fc Henningl l2004t Uohansen et al.ll2007l ). Models 
suggest grain growth to large particle sizes with a popu- 
lation of small grains remaining clo se to the surface an d 
larger grains deeper in the disc (e.g. iDominik et ai1 l2007). 
However, many details of process remain uncertain. For ex- 
ample, it is not clear how dust grains overcome the "me- 
tre size barrier" without being accreted o nto the star (e.g. 
IWeidenschillind 119771 iBrauer et all 120081) or d estroyed by 
high speed collisions (e.g. I Jones et"aT1 Il996h iBlum etlrtl 
2000), or how Kelvin-Helmholtz in stabilities prevent th e 
dust sublayer from fragmenting (e.g. Uohansen et al.| [2006) . 



Detection of the large building blocks of planets (par- 
ticle sizes > Im) in the disc midplane is far beyond cur- 
rent observational capacities, but we can expect to de- 
tect signatures of their formation. The stratified struc- 
ture resulting from grain growth and settling has di- 
rect consequences on disc observables like their spectral 
ener gy distributions (SEDs) and scattered light images 
(e.g. iDullemond fc Dominikl 12004 hereafter D04). A va- 
riety of observational techniques have been used to ob- 
tain insights int o the disc properties and their dust con- 
tent: SEDs (e.g.lD'Alessio et all 1200 ll). scattered light im- 
ages (IWatson et al.ll2007lL thermal em ission maps in the 
millimetre regi me (IDutrev et al.ll2007l). m id-infrared spec- 
troscopy (e.g. iKessler-Silacci et al.l |2006[ ). However, each 
technique only provides a limited view of a disc. SEP anal- 
ysis leads to multip le ambiguities (e.g.. lThamm et ailll994t 
IChiang et afl ! 200D since the lack of spatial resolution pre- 
cludes from solving the degeneracies between model pa- 
rameters like geometry and opacity. Spatially resolved ob- 
servations in a single spectral band (scattered light in the 
optical or near-infrared regime or thermal emission in the 
millimetre domain) also gives incomplete information about 
the disc. As the dust opacity is a steep function of wave- 
length and because high temperature gradients are present 
within the disc, a single-band observation gives insight to 
only a limited region of the disc. Thus, scattered light im- 
ages probe the surface layers of these optically thick discs 
at large radii whereas, millimetre observations are mainly 
sensitive to the bulk of the disc mass closer to the midplane. 
To precisely study fine physical processes like dust evolution 
and stratification, it is necessary to combine the aforemen- 
tioned methods in a multi- wavelength, multiple technique 
observational and modelling approach. 

In this paper, we study the circumstellar environment 
of the T Tauri star IM Lupi. The disc surrounding IM Lupi 
(Schwartz 82, HBC 605, IRAS 15528-3747) was first de- 
tected in scattered light with the WFPC2 instrument on 
board the Hubble Space Telescope (HST) as part of a 
T Tauri star imaging survey (Stapelfeldt et al. 2008, in 
prep.). IM Lupi is an M0 T Tauri star located within the 
Lupus star forming clouds. It is one of four young stellar ob- 
jects in the small 13 CO(1- 0) Lupus 2 core near th e extreme 
T Tauri star RU Lupi (|Tachihara et all II996D . Despite 
the low accretion-related spectroscopic activ ity of IM Lupi 
(|Reipurth et al.lll996l : IWichmann et alj l999). Ion ger wave- 
length observations reveal ample evidence for circumstellar 
material in the system wi th strong mm continuum emission 
dNuernberger et al.lll997l) . Single dish 12 CO and 13 CO line 
observations indicate that the dis c is gas rich and are consis - 
tent with a rotating disc model (|van Kempen et aLll2007t ). 
The 3.3 mm cont inuum emission from the disc was spa- 
tially resolved by lLommen et all (|2007f ). Preliminary mod- 
els bv lPadgett et al.f ( 20061 ) showed that the infrared excess 
is well-reproduced by a disc model. 

We have built a rich data set of observations of the disc: 
a well sampled SED, HST multiple wavelength scattered 
light images, Spitzer near- and mid-infrared spectroscopy 
and SMA millimetre emission maps. These observations 
provide different, complementary detailed views of the disc 
structure and dust properties. We investigate whether all of 
these observations can be interpreted in the framework of 
a single model from optical to millimetre wavelengths, and 
whether we can derive quantitative conclusions regarding 
the evolutionary stage of the disc. In section^ we describe 
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Fig. 1. The IM Lupi circumstellar disc observed in scattered light by HST. Left and middle panels: respectively 
F606W and F814 WFPC2 PCI PSF-subtracted images. The dark diagonal feature running through the star is an 
artifact of charge bleeding along the CCD detector columns. Right panel: NICMOS PSF-subtracted coronagraphy at 
1.6 /im. The central circle represents the 0.3" radius coronagraphic obscuration. All images are shown in square root 
stretch, with North up and East to the left. The field of view is 5x5"(dashed boxed in Fig. [2]). 



the observations and data reduction procedure. In section[31 
we first draw constraints on the disc and dust properties 
from the various observations and then, in section [H we 
present simultaneous modelling of these observations. A de- 
tailed study of the dust properties is presented in section [5] 
In section [SI we discuss the implication of the results, and 
section [7] contains the concluding remarks. 

2. Observations 

2.1. Scattered light images 

2.1.1. HST/WFPC2 observations and data processing 

HST imaging observations of the protoplanetary disc sur- 
rounding IM Lupi were obtained as part of an WFPC2 
snapshot survey of nearby T Tauri stars (HST Cycle 7 
GO/7387 program, PI: Stapelfeldt), on 1999 February 
18, using the HST Planetary Camera 1 (spatial scale of 
45. 5mas. pixel -1 ). The data consist of short and long ex- 
posures through the F606W and F814W filters (F606W: 
8 s and 100 s; F814W: 7s and 80 s) which roughly cor- 
respond to Johnson R and I. The spatial resolution in 
both filter bands F606W and F814W is 0.06"and 0.08", 
respectively. In the long exposures, the star is heavily sat- 
urated in an attempt to reveal low surface brightness cir- 
cumstellar nebulosity. After standard WFPC2 data reduc- 
tion, the stellar point spread function (PSF) was compared 
to our large database of WFPC2 PSFs, and a suitable 
match, HD 181204, was found in field position, colour, and 
exposure level. The sub-pixel registration, normalisation, 
and subtraction of the stellar PSF w as performed follow- 
ing the method of iKrist et alj (|1997[ ). Figure [1] shows the 
long F606W and F814W long exposures after PSF subtrac- 
tion. Obvious PSF artifacts remain in the subtracted im- 
age in the form of the saturation column bleed and residual 
diffraction spikes. 



2.1.2. HST/NICMOS observations and data processing 

Complementary near-IR observations were carried out on 
2005 July 11 with the Near Infrared Camera and Multi- 
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Fig. 3. Spitzer/IRS spectrum of IM Lupi. The large fea- 
tures of amorphous silicates are seen at 10 and 18 /im, as 
well as crystalline features at 9.3 /im and around 27/Ltm. 
A blend around 24 /im is tentatively detected. The spec- 
trum has been slightly smoothed to reduced the noise. The 
jump at 20 fim is an instrumental artifact resulting from 
the transition between the Long-Low and Long-High mod- 
ules of IRS. The inset shows the fit of the 10 /im continuum 
subtracted silicates feature. The thick grey line represents 
the observed spectrum (corresponding to the black line in 
the main panel) and the dot-dashed black line is the fitted 
synthetic spectrum (see section IBTTj) . 



Object Spectrometer (NICMOS) as part of the HST Cycle 
13 GO/10177 program (PI: Schneider). High contrast im- 
ages of IM Lupi's circumstellar disc were obtained using 
NICMOS camera 2's coronagraphic imaging mode (spatial 
scale 75.8 mas. pixel -1 ) following a strategy nearly identi- 
cal to that used to ima ge the debris discs around HD 32297 
(|Schneider et al J 12005] ). HD 181327 (jSchneider et alj|200tl ) 
but at a wavelength of 1.6 microns (F160W filter A c g = 
1.604 /urn, FWHM = 0.401 /im, spatial resolution of 0.16"). 
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To mitigate the potentially degrading effects of lo- 
calised, but rotationally invariant, optical artifacts that ap- 
pear in NICMOS coronagraphic images, IM Lupi was ob- 
served at two celestial orientation angles offset by 29.86°. 
IM Lupi was autonomously acquired and optimally posi- 
tioned behind the 0.3" radius coronagraphic obscuration 
following 0.241 s target acquisition imaging (also in the 
F160W filter) at each of the two spacecraft (and field) orien- 
tations. Following their respective target acquisitions, two 
sets of F160W coronagraphic observations (each yielding 
704 s of integration time after median combination of count- 
rate images derived from three STEP32 multiaccum expo- 
sures) were obtained in a single spacecraft orbit to minimise 
instrumentally induced temporal variability in image struc- 
ture on multi-orbit timescales. Direct (non-coronagraphic) 
images of IM Lupi (10.15s), were also obtained after slew- 
ing the target out of the coronagraphic obscuration. 

A set of similarly observed, high SNR, coronagraphic 
images of bright, isolated, non-disc bearing stanQ observed 
in GO/10177, served as point spread function (PSF) tem- 
plates to create PSF subtracted images of the IM Lupi 
disc. This process reduced the stellar light below the en- 
hanced contrast levels provided by the coronagraph alone. 
The PSFs of two of the template stars with J-H and H- 
K colours similar to IM Lupi, GJ 3653 (M0.5V) and GL 
517 (K5V), were well matched in structure to IM Lupi stel- 
lar PSF impressed upon the disc images. The deeply ex- 
posed F160W images of the PSF templates (total integra- 
tion times of 1352 s and 1224s, respectively, for GJ 3653 and 
GJ 517) were flux density renormalised to that of IM Lupi 
(based on target acquisition and direct imaging of all three 
stars) and used to construct four PSF subtracted images of 
the IM Lupi circumstellar disc (two from each orientation 
using both PSF templates). Residual optical artifacts aris- 
ing from the HST secondary mirror support structure were 
masked in the individual PSF-subtracted images. All four 
images were then median combined after re-orientation to 
a common celestial frame to produce the then photomet- 
rically calibrated image shown in Fig. [T] (right panel) and 
used in the analysis discussed in this paper. For further de- 
tails of the technique and process w e refer the reader to the 
afor ementioned disc i magin g papers: [Schneider et al.l ((2005) 
and lSchneider et all (|2006f) . 

2.1.3. Disc Morphology, Brightness &l Scattering Fraction 

PSF subtracted images reveal, at all wavelengths, a 
compact nebula adjacent to the SW of the star. 
Morphologically, the nebula is a broad, gentle, symmet- 
rical arc nearly 4" in size (celestial PA = 143°± 5°). No 
clear differences are found between the nebulosity at 0.606 
and 0.814 /im despite the presence of Ha and other emis- 
sion lines in the F606W filter. In the more sensitive 1.6 
micron image, the "nebulosity" is (like in the WFPC2 im- 
ages) brightest to the SW, but is clearly seen circumscribing 
the star. The scattered light to the NE is not apparent in 
the PC2 images partly due to the instrumental sensitivity 
limitations in PC2 PSF subtraction compared to NICMOS 
coronagraphy with PSF subtraction. 
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1 As determined from earlier (HST Cycle 7) NICMOS coro- 
nagraphic observations, specifically from GO/7226 for the PSF 
template stars discussed here. 



RA offset (arcsec) 

Fig. 4. SMA 1.30 mm aperture synthesis image of IM Lupi. 
Contours begin at 6 mJy/beam and step in factors of v2 in 
intensity. The black ellipse shows the CLEAN half-intensity 
beam size. The orientation is the same as in Fig. [I] The zero 
position is RA = 15:56:09.203, DEC = -37:56:06.40 Q2000). 



The IM Lupi images bear a striking morphological re- 
semblance to the protoplanetary di sc around the T Tauri 
star GM Aur ([Schneider et al.ll2003l ) and indicate the pres- 
ence of a circumstellar disc inclined in the range 40-60° 
towards the line of sight. Because of the similar appearance 
of the nebula in F606W, F814W and F160W, with a strong 
front/back asymmetry, we conclude that it is dominated 
by scattered light from the star. To the SW of the ellipti- 
cally shaped circumstellar nebulosity, an arc-like dark band 
(presumably the higher opacity disc midplane centered at 
r«1.5"along the disc morphological minor axis), bifurcates 
an isophotally concentric lower surface brightness scatter- 
ing region (presumably the " lower" scattering surface of the 
back side of the disc) extending to r«2.5"on the morpho- 
logical minor axis of the disc. This feature is most obvious 
in the F814W and F160W images. The dark lane between 
the upper and lower scattered light nebulae is best seen 
~ 1.4" southwest of the star. This coincides with the prob- 
able forward scattering direction given the purported disc 
inclination (see Fig. [2]). 

The total brightness of IM Lupi in the saturated 
WFPC2 i mages was measured using the method of 
iGillilandl (| 1994ft . The results are 0.088 ± 0.0044 Jy 
(F606W mag =11.44 ± 0.05) and 0.150 ± 0.0075 Jy 
(F814W mag =10.52 ± 0.05) in the F606W and F814W filters 
respectively. The fraction of light scattered by the disc (rel- 
ative to the stellar flux at this wavelength) is 0.0186 in the 
F814W image, measured in the region beyond 0.2"radius 
from the star. 

From the PSF-subtracted coronagraphic image, the 
measured F160W-band disc flux density at 0.3" > r > 
4.5" is 11.8mJy (uncertainty w 4%). This flux density 
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Fig. 5. Circularly averaged 1.3 mm (upper panel) and 
3.3 mm (lower panel) visibility data (crosses) compared to 
the models (lines) with different surface density exponents. 
The best fit model is shown in the full red line and corre- 
sponds to a = — 1. The dashed and dotted lines correspond 
to the best models with a — —2 and 0, where the inclina- 
tion is enforced to be 50° (see section FO)) . 

excludes the small near-central area, shown in black in 
Fig. [TJ which is insufficiently sampled due to masked op- 
tical artifacts. From the direct and target acquisition im- 
ages we measure the F160W-band stellar flux density as 
0.616 ± 0.012 Jy (F 160W may =8.11 ± 0.02) via TinyTinfl 
(|Krist fc Hooklll997t) model PSF fitting and aperture pho- 
tometry. Thus, the fraction of 1.6 /im starlight scattered 
by the disc (beyond 0.3" from the star) is 0.019 ± 0.001, 
consistent with the F814W value. 

In the F160W image, a large scale fainter halo is sensi- 
tively detected to a distance of « 4'.'4 from the central star, 
along the morphological major axis (see Fig. [5]). This large 
scale nebula may indicate the presence of a tentative enve- 
lope surrounding the circumstellar disc. Because the dark 
lane and counter nebula of the disc, which are seen through 
the potential envelope, are detected, this envelope must be 
optically thin at optical wavelengths. It is not detected in 
the F606W and F814W images. In the following, we focus 
on the disc properties and do not try to reproduce the halo 
in our modelling. 

2.2. Spitzer IRS and MIPSSED spectra 

IM Lupi was observed with the IRS spectrograph installed 
onboard the Spitzer Space Telescope as part of the "Cores to 
Discs" ( c2d) legacy progra m (AOR: 0005644800, PI: Evans, 
iKessler-Silacci et al.ll2006h . The observations took place on 
2004 August 30 using the four proposed modules (Short- 
Low, Long-Low, Short-High and Long- High), correspond- 
ing to a wavelength coverage of 5.2-38.0 /im with a spectral 
resolution R between 60-127 for the two "low" modules 
and R ~ 600 for the two "high" modules. The data re- 

2 http:/ /www. stsci.edu/software/tinytim/tinytim. html 



Table 1. Details of the Spitzer /IRS observations. SL, SH, 
LL, and LH refer to Short-Low, Short-High, Long-Low, 
Long-High respectively. R is the spectral resolution and 
SNR the signal to noise ratio. 



Module 


R 


Integration time 


SNR 


Pointing 






(tint X n^cc x n CX p ) 




offset (") 


SL 


60-127 


14 x 1 x 2 


20 


-0.8 - 0.2 


SH 


~ 600 


31 x 2 x 2 


71 


0.7 - 1.3 


LL 


60-127 


14 x 1 x 2 


35 


-1.4 - -0.3 


LH 


~ 600 


60 x 1 x 2 


66 


3.0 - 1.2 



duction was perfor med using the c2d legacy team pipeline 
(jLahuis et al.ll2006l ) with the S13 pre-reduced (BCD) data, 
using the aperture extraction methocQ. Pointing errors of 
the telescope can produce offsets between the different mod- 
ules which were corrected for by the reduction pipeline. 
Table Q] summarises the details of the observations while 
the spectrum is presented in Fig. [3] 

The 52-97 /urn MIPSSED data (R ~ 15-25) were taken 
on 2006 March 31 (Program ID 1098, PI: Kessler-Silacci). 
The basic calibrated data (BCDs) were coadded using the 
Mosaicking and Point source Extractor (MOPEX) soft- 
ware. Because the flux from IM Lupi is weak, the automated 
MOPEX extraction failed to extract the flux. MOPEX in- 
deed focused on the edge of the detector, where there are 
a lot of hot pixels, rather than on the signal from the 
source. The MIPSSED spectrum was therefore extracted 
with IRAF since it allows the user to set the centre of 
the aperture (around columns 12-16) and performs an opti- 
mised extraction. The resulting spectrum remains too noisy 
to detect any spectral features, such as the 70 /im cystalline 
emission band. In order to increase the signal-to-noise and 
get a better estimate of the mid-IR slope of the SED, we de- 
cided to bin the spectrum to calculate 3 photometric mea- 
surements at ~ 60, 75 and 90 /im (see Tabic [2]). 

The IRS spectrum of IM Lupi clearly shows silicate 
emission features typical of Clas s II objects (lHanner et al.l 
119981 : IKessler-Silacci et al J 12006). In addition to the amor- 
phous features at 10 /im (Si-0 stretching mode) and 18 /im 
(O-Si-0 bending mode), other features are visible in the 
spectrum that we attribute to crystalline forsterite (Mg- 
rich end member of the ol ivine group) and crystalline en- 
statite ('Mg-rich pyroxene. iKoike et al.| [2000: Mol ster et al.l 
2002). The 9.3 /im feature, attributed to enstatite, is clearly 
visible, as is the crystalline complex at around 27 /im. This 
latter complex is likely a blend of a forsterite feature at 
27.6 /im plus an enstatite feature at 28.2 /im (Figure [3]). A 
complex around 24 /im may also be tentatively detected at 
a signal-to-noise between 1 and 3. This complex is possi- 
bly a blend of enstatite features at 23.8 and 24.5 /tm and 
of a forsterite feature at 24 /im. The forsterite feature at 
33.6 /im is not observed. 

The shape and strength of the amorphous 10 /im fea- 
ture are related to the mean size of the emitting grains. 
These have been used as observational diagnostics for 

3 The extraction was done following two different methods: 
full aperture extraction and PSF extraction. PSF extraction is 
less sensitive to bad data samples but for some modules the 
estimated PSF is subpixel-sized. As a consequence, the PSF ex- 
traction becomes unstable. In this paper,we adopt the spectrum 
obtained with the full aperture extraction method because of its 
better stability. 
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Fig. 6. Qualitative analysis of the SED of IM Lupi. The 
blue dotted line presents a SED calculated with submicron 
particles (a max = 3 /mi). Silicates emission bands are well 
reproduced but the millimetre fluxes are too low, even with 
a disc mass of Mdi sc = 0.1 M©. A model with millimetre 
grains (a max = 3 mm, green dashed line) performs much 
better at long wavelengths but the silicate features then 
disappear. A model including millimetre grains (a max = 
3 mm) close to the midplane and submicron grains in the 
upper layers (i.e. with a stratified structure as described 
in section f3.6p . can reproduce both the silicates bands and 
the millimetre fluxes (red full line). The thin black dot-dash 
line represents the stellar photosphere. 



details and full analysis. The extended configuration of 
the SMA, with baselines up to 180 m, provided an excel- 
lent sampling of the wv-plane. The corresponding beam is 
1.84" x 1.25", with a position angle of 0.2°. 

Figure Q] presents the CLEAN map of IM Lupi's disc. 
The disc is clearly resolved with FWHMs of 2.33 x 1.72" 
and a position angle of 167°. Taking into account the beam 
elongation, this PA is in agreement with the PA of 143° 
derived from scattered light images. The emission extends 
to w 2"from the star at the 3 sigma level of 6.28mJy/beam. 

The total integrated 1.3mm flux is 175.8±4mJy. 
Together with the 3 .3 mm flux of 13mJy measured by 
lLommen et al.l (|2007| ). this leads to an estimate of the mil- 
limetre spectral index a mm = dlog(F l/ )/dlog(^) = 2.80 ± 
0.25. 

Because of the relatively large noise of individual visi- 
bilities in the wu-plane, a direct model fitting to these is not 
very informative, and some form of averaging is necessary. 
Figure [5] shows the amplitude averaged over ww-distance in 
concentric circular annuli. We compare this amplitude with 
the amplitude averaged in elliptical annuli of eccentricity 
e = 0.65 (corresponding to a axisymmetric structure ob- 
served at an inclination of 50°). The latter is aimed at pro- 
viding a more correct repre sentation of the data, taking into 
account projection effects (|Lav et al.lll997l ). The difference 
in amplitude for the two averaging methods is smaller than 
the rms levels of the observations at all baselines. Thus, 
the deprojection of the visibilities is not necessary in this 
study. In the following, we adopt the amplitude averaged 
in concentric circular annuli and use the same procedure to 
calculate the visibilities of the different models. 



grain growth in discs (e .g. Ivan Boekel et all 120051 and 
iKessler-Silacci et"aT1 [20 06). We measured the shape and 
strength of the IM Lupi 10 Mm feature f ollowi ng the ap- 
proach described in IKessler-Silacci et al.1 (|2006h . The con- 
tinuum normalisation is obtained using their first method. 
The fluxes at the reference wavelengths, S11.3 and Sg.%, 
are then calculated by integrating over a wavelength range 
of 0.1 /mi centered on A = 9.8 /mi and 11.3 /mi. Finally, 
the feature strength is estimated by calculating the mean 
peak flux, Speak, of the normalised spectrum. We obtain 
a S11.3/S9.S ratio of 1.00 and a 5p ca k of 1.60, which lo- 
cates IM Lupi in the bulk o f the class II objects observed by 
IKessler-Silacci et al.l ((2006) with a flat 10 /mi feature, con- 
sistent with the presence of micron-sized grains. Overall, 
these features indicate a significant level of processing of 
the silicate grains in the disc compared to those observed 
in the interstellar medium, which instead are found to be 
submicro n-sized and largely a morphous (<1 % of crystalline 
silicates, iKemper et al.ll2004l ). 

2.3. Spatially resolved SMA observations of 1.3 mm 
continuum emission 

The 1.3 mm (230.538 GHz) thermal dust emission from 
IM Lupi was observed with the Submillimeter ArrajQ 
(SMA). See Panic et al (in prep) for the observational 

4 The Submillimeter Array is a joint project between the 
Smithsonian Astrophysical Observatory and the Academia 
Sinica Institute of Astronomy and Astrophysics and is funded 
by the Smithsonian Institution and the Academia Sinica. 



2.4. Spectral energy distribution 

IM Lupi is an irregular variable (|Batalha et alj [l998) 
and has been obse rved to flare dramatically in U band 
(jGahm et al.l Il993h . To limit the effect of variability in 
ou r analysis, we adopt the optical photometry presented 
bv lPadgett et al.l (|2006f ). which is contemporaneous to the 
mid-IR measurements. Mid- and far-infrared flux densi- 
ties for IM Lupi w ere measured by the c2d legacy project 
(jEvans et al.l l2003). A detailed description of the reduction 
and source extraction procedures use d for the IRAC and 
MIPS measurements can be found in lHarvev et alj (S2006. 
I2007D and in the c2d delivery document ((Evans et al.ll2007l ). 
All measurements are quoted in Table O 

3. Simple estimation of some of the star and disc 
parameters 

Each of the previously presented observations gives, by it- 
self, strong insights into the properties of IM Lupi's disc and 
its dust content. Our goals in the following sections are to 
exploit this complete data set in order to draw complemen- 
tary constraints and obtain a finer understanding of the 
circumstellar environment of IM Lupi and its evolutionary 
state. We aim at building a as coherent a picture as possi- 
ble, by analysing the different observations simultaneously, 
in the framework of a single model. 

Even with simple assumptions to describe the disc struc- 
ture and dust properties, exploring the whole parameter 
space is far beyond current modelling capabilities. Instead, 
the analysis was performed in 2 steps: 
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Table 2. Photometr ic measureme nt s of IM Lupi. 

References: ( 1) = iPadgett et ail (|2006f ). (2) = 

lLommen et al. 



(120071) 



A (fj,m) 


Flux (Jy) 


Date 


Ref. 


0.545 


0.0655 ± 0.0007 


06/2004 


(1) 


0.638 


0.120 ± 0.0012 


06/2004 


(1) 


0.797 


0.216 ± 0.0022 


06/2004 


(1) 


1.22 


0.483 ± 0.0048 


03/06/2000 


2MASS 


1.63 


0.591 ± 0.0059 


03/06/2000 


2MASS 


2.2 


0.511 ± 0.0051 


03/06/2000 


2MASS 


3.6 


0.324 ± 0.0184 


12/08/2004 


Spitzer/IRAC 


4.5 


0.220 ± 0.0178 


12/08/2004 


Spitzer/IRAC 


5.8 


0.313 ± 0.0156 


12/08/2004 


Spitzer/IRAC 


8.0 


0.370 ± 0.0223 


12/08/2004 


Spitzer/IRAC 


24 


0.765 ± 0.0708 


02/08/2004 


Spitzer/MIPS 


61.1 


1.42 ± 0.22 


31/03/2006 


Spitzer/MIPSSED 


70 


1.581 ± 0.127 


02/08/2004 


Spitzer/MIPS 


74.8 


1.48 ± 0.37 


02/08/2004 


Spitzer/MIPSSED 


89.3 


1.26 ± 0.51 


02/08/2004 


Spitzer/MIPSSED 


1300 


0.1758± 0.0351 


21/05/2006 


SMA 


3276 


0.013 ± 0.0026 


08/2005 


(2) 



We extract, whenever possible, parameters "directly" 
from observations or via simple modelling: disc size from 
scattered light images, disc dust mass and maximum 
grain size from the millimetre spectral index and dust 
composition from mid-infrared spectroscopy (this sec- 
tion). We also fit the SED alone, by a manual explo- 
ration of the parameter space, to test the need for a 
spatial differentiation of dust grains within the disc. The 
goal of this preliminary analysis is to determine which 
parameters can be kept fixed, allowing us to significantly 
reduce the dimensionality of parameter space to be ex- 
plored, and which parameters need to be varied in our 
modelling. 

We then systematically explore these remaining param- 
eters, which cannot be easily extracted from observation 
and/or can be correlated with each other. We calculate 
a full grid of models and perform a simultaneous fit to 
all observations (section . 



3.1. Model description 

Synthetic images and spectral energy distributions are com- 
puted using MCFOST, a 3D continuum radi ative trans- 
fer c ode based on the Monte Carlo method (|Pinte et al.l 
2006). MCFOST includes multiple scattering with a com- 
plete treatment of polarization (using the Stokes formal- 
ism), passive dust heating assuming radiative equilibrium, 
and continuum thermal re-emission. In short, the code first 
computes the temperature structure assuming that the dust 
is in radiative equilibrium with the local radiation field. 
This is done by propagating photon packets, originally 
emitted by the central star, through a combination of scat- 
tering (following Mie theory), absorption and reemission 
events until they exit the c omputation grid. Th i s step uses 
the a l gorith ms described in lBiorkman fc Wood! (|2001[ ) and 
iLucvl (119990 with a total number of 1 280 000 photon pack- 
ets. The SEDs are then computed by emitting and propa- 
gating the proper amount of stellar and disc thermal pho- 
ton packets, ensuring that at least 12 800 packets are con- 
tributing in each wavelength bin. All of these packets are 
allowed to scatter in the disc as often as needed. We use 102 



wavelength bins, globally distributed in a logarithmic scale 
between 0.3 and 3 500/im but with an increased resolution 
(linear scale with bins of 0.5 /im) between 5 and 40 /im, to 
properly sample the regime observed with the IRS spec- 
trograph. About 4 millions packets were used to compute 
each of the scattered light images and thermal emission 
maps The maximum wavelength of HST observations be- 
ing 1.6 /im, we assume that the disc thermal emission is 
negligible in the scattered light images. For the millimetre 
maps, both the disc thermal emission and star contribution 
are considered. 

We consider an axisymmetric flared den- 
sity structure with a Gaussian vertical profile 
p(r, z) = po(r) exp(— z 2 /2 h(r) 2 ) valid for a vertically 
isothermal, hydrostatic, non self-gravitating disc. We 
use power-law distributions for the surface density 
£(r) = E ( r / r o) a an d the scale height h(r) = h {r/r^Y 
where r is the radial coordinate in the equatorial plane and 
ho is the scale height at the radius r = 100 AU. The disc 
extends from an inner cylindrical radius r- m to an outer 
limit r out . 

Dust grains are defined as homogeneous and spherical 
particles (Mie theory) with sizes distributed according to 
the power-law dn(a) oc a p da, with a min and a max the mini- 
mum and maximum sizes of grains. Extinction and scatter- 
ing opacities, scattering phase functions, and Mueller ma- 
trices are calculated using Mie theory. We adopt a power- 
law index p = —3.5, which is generally used to reproduce 
interstellar extinction curves. We choose a minimum grain 
size a m ; n , small enough that its exact value has no effect 
(we fix it to 0.03 /im). The maximum grain size a max is 
considered as a free parameter. 



3.2. Star properties 



iKrautterl (|l992f ) discussed the available distance estimates 
for the Lupus star forming regions, ranging from 130 up 
to 300 pc, and concluded that the most likely distance is 
between the limits 13 and 170 pc. The most recent de- 
terminations are from Hug hes et al.l (Il993|). who propose d 
a distance of 140±20pc and from I Wichmann et all (|l998f ). 
who concluded a distance of 190±27pc from Hipparcos par- 
allax. In the following, we adopt this value of 190 pc and 
discuss the effect of distance uncertainties on the parame- 
ters derived from our modelling in section 14.41 

IM Lupi has been class ified as an MO star from pho- 
tometric measurements by iHughes et al.l |l994). The es- 
timated age for thi s system ranges from 10 5 yr to 10 7 yr 
(jHughes et al.lll99~4l ). depending on the pre-main sequence 
evolutionary track used. In the following modelling sec- 
tions, we adopt a NextGen stellar atmosphere model 
(|Allard et al.lll997h with 3 900K and log(g) = 3.5, a stellar 
radius of 3R Q (L = 1.9 L Q ) and an Ay of 0.5 mag, which 
matches the optical pho tometry of IM Lupi. Following the 
star evolution models of iBaraffe et alJ (|l998l ). this roughly 
corresponds to a 1 M Q star with an age of 10 6 yr. 

Hipparcos observatio ns suggested a compan ion with a 
separation of only tf.'A (jWichmann et all Il998h. How ever, 
neither near-IR speckle observations ( Ghez et aLlll997| ) nor 
the current HST observations confirm the presence of a 
close companion. The extended nebulosity may have been 
marginally detected by Hipparcos and incorrectly inter- 
preted in terms of a binary. 
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3.3. Dust composition 

The strength and shape of the 10 and 20 /im emission sil- 
icate features give us some information on the dust grain 
composition. We find that micron s ized amorphous olivine 
(MgFeSiC>4, refractive indices from iDorschner et al.lll995f ) 
dust grains reproduce in a qualitative way the emission 
bands. In the following general modelling approach, we do 
not try to match perfectly the shape of the features, but in- 
stead their strength, and we fix the composition of dust to 
100 % amorphous olivine. A detailed analysis of the silicates 
features is presented in section [57Tl The results of this anal- 
ysis are consistent with the composition we adopted here. 

3.4. Disc outer radius 

Scattered light images reveal the presence of a disc with a 
diameter of approximately ~ 4". The fact that we see the 
dark lane and the corresponding counter-nebula, produced 
by scattered light from the backside of the disc, indicates 
that the disc, if it remains optically thick at optical wave- 
lengths, cannot extend much further away than 400 AU (as- 
suming a distance of 190 pc). In the following we fix the 
outer radius in our models to 400 AU. 

3.5. Disc mass 

Basically, the millimetre fluxes give access to the prod- 
uct of the disc dust mass with the dust absorption 
opacity. Assuming am orphous olivine spherical grains 
(jDorschner et all Il995f) . we obtain a maximum absorp- 
tion opacity at 1.3 mm of 2 cm 2 per gram of dust, cor- 
responding to a grain size distribution with a maxi- 
mum grain size close to 3 mm . This value is also the 
value adopted by[B cckwit h et al.l (|1990f ) and more recently 
lAndrews &: Williams! (|2007fh With this opacity, the mil- 
limetre visibilities are reproduced by a dust disc mass of 
10~ 3 M Q . This mass is estimated using the temperature 
structure from the solution of the radiative equilibrium. 
Grain size distributions with a significantly smaller or larger 
maximum grain size have a lower opacity at millimetre 
wavelengths and would imply larger disc dust masses to 
account for the observed fluxes. With the assumption that 
the gas to dust ratio is 100, as in molecular clouds, the min- 
imum total disc mass is 0.1 M Q . With a star mass close to 
IMq, the disc mass cannot be much larger than 0.1 M Q . 
We then consider it also as a maximum value and fix the 
disc mass to this value. It should be noted that millimetre 
dust opacities remain poorly known. The disc mass should 
thus be understood as an estimate with an uncertainty of 
a factor of a few. 

3.6. Assessment of the need for grain size stratification 

The SED of IM Lupi presents some clear features that we 
can use to guide our analysis on the dust properties. Our 
modelling confirms that the disc is almost entirely opti- 
cally thin in the millimetre regime. We can therefore use 
the millimetre spectral index of 2.80 ±0.25 to constrain the 
grain size and, thus the corresponding millimetre opacity 
slope is (3 mm — 0.8 ± 0.25. This value is significantly lower 
than the expected value between 1.5 and 2 for interstel- 
lar medium (ISM) grains (< 1/xm), which indicates that 
larger grains, of the order of 1mm in size, are present in 



the disc. However, the IRS spectrum clearly shows strong 
silicate emission bands that cannot be reproduced by a sim- 
ple power-law grain size distribution extending to millime- 
tre sizes, even if it contains a high fraction of micron-sized 
grains. The silicate emission requires that the maximum 
grain size do not exceed a few microns in warm (close to 
the star), optically thin regions of the disc at mid-IR wave- 
lengths (surface). In this section, we explore whether a dust 
population that is perfectly mixed with the gas, and hence 
uniform in the disc, can simultaneously reproduce both of 
these features, or whether a stratified structure with spa- 
tially separated dust populations of small and large grains 
is required to solve this apparent contradiction. 

Assuming that the dust and gas are perfectly mixed, 
the silicate emission features are best reproduced with a 
maximum grain size a max < 3 /im, essentially because larger 
grains are featureless at mid-IR wavelengths. Fig. [S] (blue 
dotted line) presents such a synthetic SED, with a max = 
3 /im. However, this model (/3 mm = 4) fails to reproduce the 
observed millimetre spectral index, as expected for grains 
in the Rayleigh regime. Furthermore, because of the low 
opacity of micron size grains in the millimetre regime, we 
also find a millimetre flux too low, by a factor of 10 relative 
to observations. 

The spectral index between 1.3 and 3.3 mm suggests in- 
stead the presence of at least millimetre-sized grains. Thus, 
a grain size distribution with a maximum grain size close 
to a max = 3 mm (Fig. [HI dashed green line) is in good 
agreement with our millimetre observations. This kind of 
grain size distribution however, does not reproduce the 
silicate emission features which almost completely disap- 
pear. Indeed, for a grain size distribution with a slope 
p = —3.5 and a maximum grain size > 100 /mi, the mid- 
IR opacities are not dominated by the sub-micron sized 
grains but rather by the grains a few microns in size 
which have featureless op acities (see for instance Fig. 5 of 
iKessler-Silacci et al.ll2006h . We can then conclude that mil- 
limetre grains are present in the disc, but that they do not 
dominate the opacity in regions producing silicate emission 
features: inner edge and disc surface at radii smaller than 
a few astronomical units. A spatial dependence of the dust 
properties, particularly grain size distribution, seems to be 
required to account for all observables. In the following, we 
assume that this spatial dependence is produced by verti- 
cal dust settling in the disc. This structure is characterised 
by small grains on the disc surface, which produce the sili- 
cate emission, and larger grains in the midplane which emit 
most of the thermal radiation in the millimetre regime. We 
discuss other possible explanations for the grain size strati- 
fication in section We describe this dust stratification by 
a simple parametric law, the reference scale height being a 
function of the grain size 

h (a) = h (a min ) ( — — J . (1) 

In the case where the dust and gas are perfectly mixed, 
£ = 0, and h {a) = /i («min) is independent of the grain 
size and is equal to the scale height of the gas disc. 

The full red line in Fig. [6] presents a SED calculated 
with a max = 3 mm and the stratification described above, 
taking £ = 0.1. This kind of model reproduces both the 
silicate emission bands and millimetre fluxes and gives sup- 
port to the presence of a stratified structure in the disc of 
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Fig. 7. SEDs of the best models. The dashed red line cor- 
responds to best model for the fit of the SED only and 
the full green line corresponds to the simultaneous fit of all 
observations. 



IM Lupi. We would like to emphasise that, at this point, 
this model is not a fit, just an educated guess as far as 
disc parameters are concerned, meant only to illustrate the 
impact o f sett ling on the SED. A simil ar effect was also 
found bv lD04l and lD'Alessio et~aT1 (|2006f) . In the following, 
we fix the maximum grain size to a max = 3 mm to match 
the spectral index and consider £ as a free parameter. 



4. Inferred disc model parameters from 
multi-technique model fitting 

Encouraged by the previously described qualitative be- 
haviour, we then moved to a quantitative analysis by ex- 
ploring a full grid of models. The goal is to find a single 
model which self-consistcntly fits all observables and esti- 
mates the parameters as well as their range of validity and 
uncertainties. Because the F606W and F814W images are 
very similar, we restrict our fitting of the scattered light im- 
ages to the F606W and F160W bands. Since they are the 
shortest and longest wavelengths: the intermediate F814W 
image is likely to be well fitted by a model that adequately 
represents these 2 wavelengths. 

4.1. Explored parameter space 

Exploration of the parameter space was performed by 
varying the geometrical properties of the disc (r- m , /?, a, 
/io(a m in)) and the degree of the dust settling £. Scattered 
light images give us a good idea of the inclination angle, 
around 50°, but we consider it as a free parameter in or- 
der to fine tune our estimation. This results in a total of 6 
independently varied parameters. The range explored for 
each parameter is summarised in Table |3] All combina- 
tions of parameters were explored, resulting in a total of 
421 200 calculated models. The ranges were chosen to sam- 
ple physically plausible values for the different parameters, 
and adapted during the course of the modelling to calculate 
all the models which may be a reasonably good represen- 
tation of data. The probability curves we obtain (Figure [9|) 
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for the different parameters show that we have sampled a 
large enough fraction of the parameter space. 



4.2. Fitting procedure 

The fitting of the SED is performed through the definition 
of a reduced \ 2 : Xsed h' om the observed and synthetic 
fluxes, and from the observational uncertainties. Because 
synthetic SEDs are sampled with a large number of packets, 
we neglect the model Monte Carlo error bars which are 
significantly smaller (< 1%). We can get a sense of their 
amplitude in the silicate feature (inset in Fig. [7]) . 

Comparing synthetic images with observations is a more 
delicate procedure due to the presence of the star that dom- 
inates the flux in the synthetic images. Observational ef- 
fects are taken into account by choosing pixel sizes similar 
to those of observations, i.e. 0.045 "/pixel at 0.6 (im and 
0.076 "/pixel at 1.6 p,m. Real scattered light images were 
obtained by manual subtraction of the PSF. Such a proce- 
dure cannot be considered for the models however, given 
the large number of models dealt with here. To make the 
problem of PSF subtraction tractable, we choose to con- 
volve our models only with the PSF core (up to 5 pixels 
from the peak) in order to reproduce its smoothing effect, 
but at the same time avoiding the superimposition of the 
convolved direct stellar light on the disc, where it is de- 
tected in the actual PSF-subtracted images. The disc ap- 
pears left-right symmetrical in both WFPC2 and NICMOS 
images. To increase the signal-to-noise ratio and to reduce 
artifacts, the observed images were symmetrised relative to 
the semi-minor axis of the disc, prior to comparison with 
models. 

A pixel-by-pixel comparison between models and obser- 
vations is very sensitive to observation artifacts. It requires 
a precise map of the observational uncertainties in the im- 
age in order to make optimal use of its information. This 
procedure is complex for high contrast PSF subtracted im- 
ages. Large deviations in small regions of the image can 
strongly bias the fitting procedure. Also, PSF subtraction 
uncertainties introduce systematic, as opposed to random, 
errors, which are not properly taken into account in a x 2 
minimisation. Instead, we adopt an image fitting performed 
by extracting geometrical observables. Our goal here is to 
reproduce the main characteristics of the images, allow- 
ing small parts of the model images to not be in perfect 
agreement with the observations, if necessary. We therefore 
adopt the following method, illustrated in Fig. [SJ which is 
based on the azimuthal intensity variations. For both wave- 
lengths, we calculate the observed intensity variations nor- 
malised to the average flux in the azimuthal angle range [- 
15°, +15°] (centered on the disc's semi-minor axis). For each 
10° sector, we extract the average flux of the pixels encom- 
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Fig. 8. Scattered light images of the best models compared to observations. Left panel: The upper row corresponds to 
the images at 0.606 /im and the lower row to those at 1.6/xm. Synthetic maps (right) were convolved by the core of the 
PSF (up to 5 pixels from the peak). Central panel: 0.6 /im azimuthal brightness profile. Right panel: 1.6 /im azimuthal 
brightness profile. In both central and right panels, the red dashed line correspond to the best fit of the scattered light 
only, the full green line to the best fit of all observations simultaneously and the blue dotted to the scattered light images 
with porous grains. The azimuthal angle is 0° in the front side of the disc, i.e. towards bottom in the first panel. 



passed between two ellipses of eccentricity e = 0.65 (cor- 
responding to a circularly-symmetric structure observed at 
an inclination of 50°) and of semi-major axis of 1.5 and 
2.1". Regions presenting diffraction artifacts are moreover 
excluded from these sectors. Uncertainties were calculated 
by repeating the same procedure on the observations from 
which we have added and subtracted the estimated error 
maps. These error maps are constructed by taking the shot 
noise from photon statistics and adding in a PSF scaled by 
the uncertainty factor in the normalisation used for PSF 
subtraction. 

The synthetic brightness profiles were extracted in the 
same way from Monte Carlo images. We do not consider any 
uncertainties for these synthetic profiles. For the 1.6 /im im- 
age we further use the positions of the dark lane and second 
nebula as geometrical observables to fit for. We compare 
the synthetic maps and the observations using reduced x 2 
based on the previously defined observable^: Xo 6 am an< ^ 

2 

A 1 . 6 fim ■ 

Model comparison with SMA and ATCA data was per- 
formed in (u, v) plane to avoid additional uncertainties 
from the image reconstruction with the CLEAN algorithm. 
Synthetic visibilities were computed by Fourier transform 
of the MCFOST emission maps. The visibilities were then 
binned in circular annuli as done for the observed visibili- 
ties, and the reduced \ 2 is computed from these binned vis- 
ibilities: Xmm- Because only 6 data points are available for 
the ATCA observations, it is not possible to define a mean- 
ingful reduced x 2 for these observations alone. Instead, we 
fit simultaneously the ATCA and SMA data by defining a 
unique reduced x 2 - To avoid an eventual contamination by 
a potential envelope and/or surrounding cloud, fitting was 
only performed for baselines larger than 40 kA which probe 
spatial scales smaller than 5 ", i.e. corresponding to the disc 
observed in scattered light. 

We compute a total Xtot defined as the sum of all the 
reduced x 2 : 

2 _ 2 _,_ 2 _j_ 2 _i_ 2 (n\ 

Xtot — XSED > Xo .6 urn ' Xl.6 /im ' Xmm - \ ) 



5 A total of number of 31 and 34 measurements were used at 
0.6 and 1.6 /im respectively. 



This Xtot allows us to do a global fitting of all observations 
simultaneously. 

4.3. Best model 

The parameters of the model with the lowest Xtot are de- 
scribed in Table |U Overall, the best fit model is in very 
good agreement with all observations, especially taking into 
account the range of wavelengths and the variety of obser- 
vations analysed in the fitting procedure. Table[5]shows the 
X 2 corresponding to best models fitting only one of the ob- 
servations and to the global best model. The global best 
models has x 2 values only slightly larger than the x 2 of the 
best models fitting each of the observations, showing that 
it offers a good representation of all observations. 

The SED of this model is represented in Fig. \7\ It is 
in very good agreement with observations from the optical 
to the millimetre regimes. In particular, it reproduces both 
the silicate bands and millimetre fluxes, and hence the mil- 
limetre spectral index. The silicates features are roughly 
reproduced but the shape is not exactly identical to the 
observed one, resulting in relatively high values of the x 2 - 
This is due to our very simple approximation on optical 
properties, with dust grains only composed of amorphous 
silicates. We recall here that a precise match of the silicate 
emission bands was not our goal at this stage of the mod- 
elling and that we were only interested in the amplitude of 
the silicate features. 

Synthetic images of the best model are compared with 
observations on Fig. [5J The general shape is well reproduced 
at both wavelengths, with the right roundness, brightness 
distribution and a dark lane in good agreement with obser- 
vations. The central and right panels show a quantitative 
comparison of the azimuthal brightness profiles. The agree- 
ment is good as shown by the values of the x 2 m Table [5l 
The model predicts profiles that are very similar at both 
wavelengths whereas the observations suggest (given the 
relatively low SNR and large PSF uncertainty) a stronger 
contrast at 0.6 /im, with a back side compatible with a non- 
detection. We discuss this point in detail in section I5T21 

Figure [5] shows the results of the fit of the 1.3 and 
3.3 mm visibilities. The agreement is also very good for the 
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baselines > 40 kA on which the fitting was performed. At 
smaller baselines, the best models are marginally below the 
data points. This may be due to a small contamination of 
the observations by the surrounding molecular cloud and/or 
envelope. Millimetre visibilities are mainly sensitive to the 
apparent spatial distribution of the dust, i.e, the surface 
density exponent a, the disc outer radius, the total mass 
and the disc inclination. With an outer radius of 400 AU de- 
rived from scattered light images and a dust disc mass based 
on the millimetre flux, our modelling allow us to constrain 
the surface density profile (Fig. [5]). The a=-l model (full 
line) provides an excellent match. Conversely, the other two 
models (a—0 and a— -2) fail to successfully reproduce the 
data, if we force the inclination to be close to 50°. Models 
with a=0 provide a good match to the observations for a 
disc closer to edge-on (inclination of « 70°), which is in- 
compatible with the scattered light images. 

Table 4. Best model parameters. The best model is the 
model with the lowest Xtot- The valid range for each pa- 
rameter is defined as the range of values which symmetri- 
cally enclose the central 68 % of the probability (see text 
for details). This is equivalent to a 1 a confidence interval. 



parameter best model valid range 



Tcff (K) 


3 900 


fixed 


-Rstar (Rq) 


3.0 


fixed 


distance (pc) 


190 


fixed 


rout (AU) 


400 


fixed 


Mdust (Mo) 


10" 3 


fixed 


dust composition 


100 % amorphous olivine, fixed 


ttmin 


0.03 /im 


fixed 




3 mm 


fixed 


<(°) 


50 


45 - 53 


surface dens, exp 


-1 


-0.84 - -1.21 


flaring exponent 


1.15 


1.13 - 1.17 


r in (AU) 


0.32 AU 


0.25 - 0.4 


ho (&min) 


10 AU 


9.5 - 10.3 


settling exponent £ 


0.05 


0.02 - 0.07 



In Fig. [7| and we present both the overall best fit 
model and best fits to single observables. This quantita- 
tively shows that the best overall fit is very similar to the 
individual fits and illustrates how limited single data sets 
arc: families of models can be found that give equally good 
fits to the observations. It is only through a combination of 
the various observations that we can solve most ambiguities 
and disentangle between models (Fig. [5]) . 

4.4. Validity range of parameters 

To determine the range of validity of the different parame- 
ters, we use a Bayesian inference method ()Press et al.lll992t 
lLav et al.|[l997tlPiiite et al.ll2007l ). This technique allows us 
to estimate the probability of occurrence of each value of a 
given parameter. The relative probability of a point of the 
parameter space (i.e. one model) is given by exp(— x 2 /2) 
where x 2 refers to the previously defined reduced x 2 of the 
corresponding model. All probabilities are normalised at 
the end of the procedure so that the sum of the probabil- 
ities of all models over the entire grid is equal to 1. This 
method analyses all the models in a statistical sense and 
does not give any specific attention to a given model, in- 
cluding the best model. 



The Bayesian method relies on a priori probabilities on 
the parameters. Here, we assume that we do not have any 
preliminary available information, and we choose uniform 
a priori probabilities which correspond to a uniform sam- 
pling of the parameters. However, in the absence of any 
data, some grid points are more likely than others: consid- 
eration on solid angles show that an inclination between 
80 and 90° (close to edge-on) is more likely than a inclina- 
tion between and 10°. Uniformly distributed disc inclina- 
tions and orientations in three dimensions correspond to a 
uniform distribution in the cosine of the inclination. Also, 
some physical quantities, like the inner radius, tend to be 
distributed logarithmically. The grid was built according to 
these distributions (Table [3]). 

Figure [9] presents the relative figure of merit estimated 
from the Bayesian inference for each of the parameters, for 
the fitting of the scattered light images, SED, and millime- 
tre visibilities. These results were obtained by marginalising 
(i.e. summing) the probabilities of all models where one pa- 
rameter is fixed successively to its different values, i.e. the 
probabilities of our 6-dimension parameter space are pro- 
jected successively on each of the dimensions (they are not 
cut through the parameter space). Potential correlations 
between parameters are entirely accounted for with this ap- 
proach, and the error bars extracted from the probability 
curves take into account the interplay between parameters. 
The resulting histograms represent the probability that a 
parameter takes a certain value, given the data and assump- 
tions of our modelling. 

Both scattered light images are similar and naturally 
lead to similar probability densities (blue dashed and red 
dotted lines). Because the error bars are smaller for the 
NICMOS image, the resulting constraints on the parame- 
ters are stronger. The disc inclination is well constrained 
in the range 45 - 55° . Fitting of both images tends to pre- 
fer flaring exponents > 1.1. Because the disc is optically 
thick at optical and near-infrared wavelengths, scattered 
light images are insensitive to the surface density. These im- 
ages probe the outer part of the disc, and due to the flared 
geometry adopted in the current work, the stellar illumina- 
tion does not depend on the inner radius. As a consequence, 
we do not get any constraint on the inner edge of the disc. 
The azimuthal brightness profile is strongly sensitive to the 
disc scale height. Both images give slightly different, but 
compatible, constraints on the scale height, with a peak of 
probability at 10 AU. The probability curve for the 0.6 /im 
image also presents secondary peaks at 5 and 15 AU. These 
peaks correspond to the two inclination bins adjacent to 
the most probable bin centered at cos(i) = 0.65. 

The analysis of the SED gives different constraints 
(green dot-dash line in Fig. [9]). As expected, the inclination 
is not well constrained. Only the models with high inclina- 
tion, i.e. models that occult the star, are excluded. As the 
disc is very massive, it remains opaque to its own radiation 
up to long wavelengths, and the SED is only slightly sen- 
sitive to the disc surface density index. The flaring index, 
describing the disc geometry and its capacity to intercept 
stellar light, is strongly constrained with most probable val- 
ues between 1.1 and 1.15. The inner radius is constrained 
between 0.15 and 0.6 AU. The probability curve for the 
settling index is of particular interest. The value 0, cor- 
responding to the case without settling, has a significantly 
lower probability than those corresponding to a relatively 
low amount of settling. The peak of probability is around 
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Fig. 9. Bayesian probabilities of the various parameters for the scattered light images at 0.6 ^m (dashed blue) and 1.6 fim 
(dotted red), the SED (dot-dash green), the millimetre visibilities (dot-dot-dash pink) and for the images, SED and mm 
visibilities simultaneously (full black line). The triangles represent the parameters of the best model. 



£ = 0.025 and is followed by a decreasing probability den- 
sity towards strongly settled discs. The probability becomes 
negligible for £ > 0.15. Indeed, high degrees of settling re- 
sult in silicate features that are too strong and mid-infrared 
fluxes that are too low. The scale height is constrained to 
be around 10 AU, in agreement with the modelling of the 
scattered light images. 

The pink dot-dot-dash line shows the constraints re- 
sulting from the fitting of the millimetre visibilities. As ex- 
pected, because they are mainly sensitive to the bulk of the 
disc mass, these observations do not give any constraints 
on the inner radius, scale height or degree of dust settling. 
Because our fitting was realised on azimuthally averaged 
visibilities, the inclination is not well constrained. Only the 
close to edge-on disc models are excluded. As previously 
mentioned however, the CLEAN map does give an inclina- 
tion consistent with the one deduced from scattered light 
images. The main constraints are obtained on i) the flaring 
index, with the most probable value between 1.1 and 1.2, 
which is in agreement with the results of the SED modelling 
and ii) on the surface density index, with a peak of proba- 
bly close to —1, and the extreme values of and —2 being 
excluded. There is a strong correlation between the sur- 
face density index and the inclination, models with a > — 1 



corresponding to inclinations > 55°. Indeed, the millimetre 
interferometers probe the apparent spatial distribution of 
the dust and an inclined disc can mimic a disc with a more 
concentrated dust distribution. 

Combined, the various observations give complemen- 
tary constraints on the model parameters. Figure [5] illus- 
trates, in a quantitative way, the complementarity between 
the results extracted from the scattered light images, the 
SED, and the millimetre visibilities. The Bayesian analysis 
gives the relative likelihood pi(M\Di) of different models 
M, given the data Di (scattered light images, SED or mil- 
limetre visibilities) that have been measured. Given the n 
different data sets, the relative likelihood of a model is then: 

n 

p(M\D 1 ...D n ) = l[p i (M\D i ), (3) 

i=l 

because the uncertainties in each data set are independent. 
In our case, this is equivalent to calculate the probability 
from a Xtot we have defined as the sum of the individual 
reduced \ 2 - The black line in Fig. [5] shows the resulting 
marginalised probabilities, corresponding to the simultane- 
ous fitting of the two scattered light images, the SED, and 
the millimetre visibilities. All parameters are constrained 
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Table 5. Reduced x 2 values for the best models. 



Model 


Xo.6 /^m 


Xl . 6 f-im 


XSED 


Xmm 


Xtot 


Best 0.6 /an 


1.16 


17.24 


83.47 


83.60 


185.47 


Best 1.6 /jm 


3.02 


0.48 


29.75 


18.64 


51.89 


Best SED 


22.83 


139.53 


9.87 


51.24 


223.48 


Best mm 


22.69 


46.13 


560.75 


0.43 


630.00 


Best all obs. 


1.56 


0.68 


14.02 


1.18 


17.44 



within narrow ranges and the corresponding probability 
curves are much sharper, i.e. the constraints on the param- 
eters are stronger (see Table 2]). This illustrates the need 
for simultaneous modelling of various kinds of observation 
to quantitatively derive disc parameters and obtain finer 
models. Table 0] gives the range of validity of the different 
parameters. For each parameter 9 with a density of prob- 
ability p(8), this range of validity is defined as the interval 
62] where 



P {9) d9 



p(9) d9 



(4) 



with 7 = 0.68. The interval \9\, 6*2] is a 68% confidence 
interval and corresponds to the 1 a interval for a Gaussian 
density of probability. 

The triangles in Fig. [9] represent the different parame- 
ters of the best model. Although this is not necessarily the 
case in Bayesian studies, it coincides here with the most 
probable model, defined as the model with the parame- 
ters which have the highest probabilities. This means that 
the best model corresponds to the global minimum of x 2 in 
the parameter space. There is no contradiction between the 
constraints from the different observations indicating that 
our model, although based on simple assumptions, provides 
an adequate description (at the precision level of present 
observations) of the disc structure and dust properties. 

4.5. Effect of distance 

The results described thus far have been determined using 
a distance of 190 pc. It should however, be noted that the 
distance of IM Lupi is relatively uncertain (see section [3~2f . 
An alternate distance value will thus affect the results of 
our modelling. Fortunately, due to the self-similarity of the 
radiative transfer equations, our results can be scaled for 
an alternate distance. If instead the distance to IM Lupi is 
A x 190 pc, where A is a real positive number, the calcula- 
tions will be mathematically identical to the ones that we 
have presented if (i) all lengths in our model (star radius, 
disc inner and outer radius rj n and r out ) are multiplied by 
A, (ii) the geometry is not modified (the scale height at 
100 AU, flaring and surface density indices, degree of dust 
settling remain identical) and (iii) the disc mass is multi- 
plied by A 2 (ensuring that the disc opacity is not modified). 



5.1. Disc mineralogy 

5.1.1. Compositional fitting of the 10 /im silicate feature 

A detailed study of the 10 /im feature provides quan- 
titative information on the composition and size of 
the grains responsib l e for the observed emission (e.g. 



Bouwman et al. 2001 



."EL 

20071 IBouv et alj 120081 and iBouwman et al.ll2008l ). To fit 

the 10 /xm fea ture, we use five dust s pecies: olivine (glassy 
Mg Fe SiC>4 , IDorschner et al.l 19951 ) , pyroxene (glassy 
MgFe SiC>6, IDorschner et al.l 11995) and silica (amorphous 



2005. Schegerer et al 



van Boekel et ail 120051 lApai et al 
20061 or more recently iMerin et al.l 



quartz (silicon dioxide) at 10K, Henning & Mutschke 1997 
for the amorph ous grains, and enstatite (jJaeger et al.l ll998) 
plus forsterite ([Servoin fe Pirioul I1973D for the crystalline 
species. Furthermore, we consider two single representative 
grain sizes (generally 0.1 and 1.5 /im, although other grain 
sizes have been tested) following previous studies. Once the 
continuum emission has been subtracted, a fit to the amor- 
phous 10 /<m feature is calculated. This is done by multi- 
plying grain opacities with relative mass fractions and with 
a single-temperature blackbody. A pseudo-x 2 minimisation 
procedure provides the best relative mass fractions for all of 
the species and grain sizes, as well as the best temperature. 

Achieving a good estimation of the dust continuum 
emission contributing to the spectrum is the most difficult 
part for this sort of analysis, because it influences the de- 
rived composition. We learned from previous studies that a 
simple local continuum of the 10 /im feature is not a good 
solution, despite the fact that it is the simplest one. We 
find that it is necessary to leave flux at the end of the 
amorphous feature, at around 13.5 /im, because amorphous 
silicate grains are still contributing to the emission at these 
wavelengths. We therefore adopt the continuum produced 
by the MCFOST code to obtain a more realistic estimate of 
the continuum shape around 10 /im. Nevertheless, the best- 
fit disc model includes amorphous silicate grain emission, 
and consequently the SED is not featureless. We generated 
another continuum by artificially removing the amorphous 
features in the MCFOST calculations. Refractive indices 
in the wavelength range 7-35 /im, which includes the sil- 
icates features, were replaced by a log — log interpolation 
with wavelength of the refractive indices from the values 
at A =7 and 35 /im. This solution, although providing a 
more physical continuum, remains imperfect. First, the in- 
terpolated refractive indices are very likely to be a strong 
simplification of what would be the refractive indices with- 
out silicate features. Second, the energy normally escaping 
through the amorphous features is redistributed over adja- 
cent wavelengths and this overestimates the flux near the 
feet of the amorphous features. A smoothing treatment is 
thus required to obtain a good estimate of the continuum. 

5.1.2. Results 



5. Detailed analysis of the dust properties 

In the previous section, we performed a global fit to all of 
the observations simultaneously. Such a fit however cannot 
account for the fine signatures we detect in the various ob- 
servations. In this section, we analyse in detail the silicate 
emission bands and scattered light images to infer addi- 
tional information on the dust properties. 



We investigated the influence of the grain size on the good- 
ness of the fit, the result of which is that only grains of 
an approximate size of 1.5 /tm are able to reproduce the 
amorphous feature (inset in Fig. [3]). A first fit with grain 
sizes of 0.1 and 1.5 /im reveals that 94% of the grains need 
to have a size of 1.5 /tm. A second run assuming grain sizes 
equal to 1.5 and 3.0 /tm indicates that 97% of the grains 
have a size of 1.5 /tm. Because 3.0 /im grains are not feature- 
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less at 10 /im, our results on the low fraction of these grains 
are robust. This confirms the qualitative result we obtained 
with the ratio Sn.s/Sg.s versus the Speak, which indicated 
micron-sized grains responsible for the 10 /iin feature. The 
result is very robust against using different continua. The 
derived temperature of these 1.5 /im grains is always around 
350 K. The best fit is obtained for a temperature of 357 K. 

In the radiative transfer modelling of the disc, the flux 
received by the observer in the 10 /im silicate band mainly 
comes from grains slightly smaller than 1 /im in size : 50 % 
of the energy is emitted by grains between 0.1 and 1.2 /im 
and the the peak of the emission originates from grains 
0.9 /im in size. This explains the slightly sharper silicate 
features in the model compared to observations (Fig. [7]). 
The average temperature of the grains responsible for the 
received emission at 10 /im is 690 K, which is significantly 
higher than the value for the fit of the silicate feature. This 
may indicate that grains around 1.5 /im, responsible of most 
of the emission in the silicate band, are located at larger 
radii and/or in deeper regions than what is assumed in 
the global modelling. However, the large difference between 
temperatures should be interpreted with care and should be 
considered indicative only. Indeed, contrary to the global 
modelling, the detailed composition fit of the silicate fea- 
ture has been performed after subtraction of the continuum 
from the MCFOST calculations (then at a temperature of 
690 K). This continuum represents about 50% of the flux 
at a wavelength of 10 /im. 

In order to reproduce the amorphous silicate feature, 
the fraction of crystalline grains is always less than 10%. 
For the best fit, the crystallinity fraction is 7%. This can 
mainly be explained by the presence of the enstatite fea- 
ture at around 9.3 /im. Long-wavelength features, due to 
enstatite and forsterite grains, also indicate the presence of 
crystalline grains in cooler regions of the disc (larger radii 
and /or deeper in th e disc) than the regions probed at 10 /im 
(jMerfn et al.ll2007l ). However, these features remains weak. 
No features are detected at wavelengths larger than 30 /im 
which also suggests a low degree of crystallinity. 

5.2. Are we observing fluffy aggregates? 

Light scattering strongly depends on the dust properties, 
with the general trends that grains significantly smaller 
that the wavelength produce isotropic scattering and that 
the larger the grains the more anisotropic the scattering 
phase function becomes. All of IM Lupi's scattered light 
images present a strong front-back asymmetry indicating 
anisotropic scattering. This implies the presence of grains 
at least comparable in size with the wavelength, in the re- 
gions of the disc probed by scattered light, i.e. the disc 
surface at large radii (> a few 10 AUs). Interestingly, the 
front-back flux ratio is higher at the shortest wavelength 
indicating that scattering could be more forward throwing 
in that case. Our modelling has lead to solutions that are 
marginally compatible with the flux level in the back side 
of the F606W image, although with a systematic overes- 
timation (Fig. [8}. In this section, we explore what kind 
of dust grains could produce a more contrasted azimuthal 
brightness profile at 0.6 /im, whilst remaining in good agree- 
ment with the azimuthal brightness profile at 1.6 /im, i.e. 
we try to find dust properties that could produce a flux 
ratio F(180°)/F(0°) of 0.05 at 0.6 /im and of 0.2 at 1.6 /im 
(180° and 0° refer to the azimuthal angles as in Fig. [HI and 



not to the scattering angles). As our best model is already 
marginally compatible with the observations, we do not try 
to perform a new fit to all observations simultaneously, in- 
stead only focusing on the effect of the dust properties on 
the flux ratios. We adopt the geometry of the previously es- 
timated best model. We do not consider dust settling here 
and take a max as a free parameter. This value should now be 
understood as the maximum grain size in the disc surface. 

The scattering anisotropy can be described, as a first 
approximation, in terms of the asymmetry parameter g\ = 
(cos 9\) where 9\ is the scattering angle at the wavelength 
A. The models we have presented so far correspond to g val- 
ues around 0.6 for both wavelengths, whereas for a better 
agreement on the flux on the backside of the disc, models 
with 5o.6/jm ~ 0.8 are needed. In Fig. [TO) we plot the re- 
gion (50.6pm ~ 0.8, (/i.6pm ~ 0.6) as a shaded area. This 
corresponds to models that would reproduce both observed 
azimuthal brightness profiles. 

The dust grains we have used so far cannot account for 
the observed scattering properties, regardless of the maxi- 
mum grain size fFig. 1101 full line). Indeed, they cannot re- 
sult in an asymmetry parameter larger than ss 0.65. Within 
the assumption that the scattering properties of the grain 
can be represented with the Mie theory, there are two ways 
to increase the scattering anisotropy: larger grain sizes and 
lower refractive indices. 

Then, we first tried to see whether we could reproduce 
the g values at both wavelengths with a different grain size 
distribution, by varying the slope of the distribution. We 
find a solution with a slope p = (Fig. [101 dot-dash line), 
which is noticeably flatter than that of interstellar medium 
grains. This corresponds to a dust population almost de- 
void of ve ry small grains (< 0.1 p). Models of dust co- 
agulation (jDullemond fc Dominikll2005t lOrmel et al"1l2007h 
show that, at least during the first stages of the process, 
the slope of the grain size distribution does not vary much 
in the regime that we are interested in. Therefore, we con- 
sider that the solution with a flat grain size distribution is 
unlikely to occur. 

The optical properties of amorphous olivine we have 
used correspond to compact grains. An attractive solution 
to increase the scattering anisotropy is to consider porous 
grains, composed in a high fraction by vacuum, which will 
result in lower refractive indices. The porosity is defined 
as the fraction of the grain volume composed of vacuum: 

*P ^vacuum/^grain — 1 ^solid/^grain, where Vvacuum, ^solid 

and V^ ra i n represent the volumes of the vacuum compo- 
nent, of the solid component (olivine silicates in this case) 
and the total volume of the grain, respectively. We calcu- 
late the optical properties of these porous grai ns assuming 
the B ruggeman effective medium mixing rule (|Bruggemanl 
Il935t ). The dashed line in Fig. [TU] represents the multi- 
wavelength scattering properties of olivine grains with a 
porosity V — 0.8. These grain complexes, with a maximum 
grain size of around 0.7 /im, mimic well the scattered prop- 
erties inferred from HST images and in particular predict a 
more forward throwing scattering at 0.6 /im than at 1.6 /im 
(Fig. [8l blue dotted lines). We consider this as an indication 
that we may be observing scattered light by highly porous 
dust grains. 

We also explored other dust compositions, in par- 
ticular the porou s mix ture of silicates and carbon of 
iMathis fc Whiffenl (|l989l model A). This model also has a 
porosity V — 0.8, and the scattering properties are very 
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Fig. 10. Asymmetry parameter as a function of maximum 
grain size for different dust compositions. The top panel is 
for A = 0.6 fim and the bottom panel for A — 1.6 fim. The 
full line corresponds to the olivine compact grains used in 
the global modelling (slope of the grain size distribution 
p = —3.5). The dashed line corresponds to the same grains, 
with a porosity of 80% (p = —3.5). The dot-dashed line 
corresponds to compact grains with p = 0. The dotted line 
corresponds to pure water ice grains. The shaded areas rep- 
resent the region that would fit the flux ratio of both the 
F606W and F160W scattered light images. A model repro- 
ducing both scattered light images must lie on the shaded 
area of both panels for the same maximum grain size. 

similar to those of our porous olivine dust grains. This 
model only predicts very faint silicate emission features and 
cannot account for the observed IRS spectrum. However, 
this reinforces our conclusion that grains with high poros- 
ity, independently of their exact composition, may provide 
a good explanation for the observed brightness profiles. 
Interestingly, scattering properties of porous grains have 
been found to be a good representation of those of fluffy 
aggregates, as lo ng as the size of the inclu sions is in the 
Rayleigh regime (jVoshchirinikov et al.l [2005) . If the inclu- 
sions are larger, then the accuracy of the effective medium 
theory becomes insufficient and more complex models, such 
as multi-layered spheres, are required to obtain a precise de- 
scription of scattering properties. But the general trend, of 
an increase of forward throwing scattering with porosity, 
remains valid. 

The presence of porous aggregates suggested by the 
scattered light images is very likely to also modify the 
SEDs, in particular the profile of the silicate bands, shift- 
ing th e peak wavelength and affecting the shap e of the 
bands (jMin et alJ200alVoshchinnikov fe Hennindl2008h. as 
well a s the dust emission in the millimetre regime ([Wright! 
Il987l ) . B oth the opaci ty at long wavelength (see for instance 
Fig. 4 in Wri ghtl[l987f ) and in the silicate bands are however 
strongly dependent on the shape of the aggregates, and not 
only o n the degree of porosity. IVos hchinnik ov fe Henningl 
(2008| find that sm all porous grains have signatures of 
large grains whereas iMin et all (|2006l ) find that large ag- 
gregates composed of small spheres have signatures of 



small grains. Different methods are used to represent the 
aggregates in both cases (multi-layered spheres and dis- 
crete dipole approximation, respectively) indicating that 
the micro-structure of the grains is very likely to be a key 
element in the opacity of aggregates. 

An alternative explanation is the presence of ice man- 
tles (H2O, CO, . . . ) around dust grains. Ices have refractive 
indices significantly smaller than rocks 1.3 for water ice 
and w 1.2 for CO2 for instance). The dust grain tempera- 
ture at the surface of the disc, for radii larger than 100 AU, 
does not exceed 70 K, allowing the formation of such ice 
mantles. The dotted line in Fig. [TOl prese nts the scattering 
prop erties of grains composed of water ice ([Irvine fc Pollack! 
Il968l ). With a maximum grain size around 0.9 /xm, these 
grains can also reproduce the scattering anisotropy at both 
wavelengths. 

Current observations do not allow distinction between 
these two solutions (porous rock grains or grains with 
ice mantles). A s sho wn in iGraham et al.l |2007) and 
iFitzgerald et al.1 ([2007) in the case of the debris disc sur- 
rounding AU Microscopii, resolved polarimetry is a pow- 
erful tool to solve this ambiguity. Only porous grains can 
produce strongly anisotropic scattering and a high level of 
polarization. 

6. Discussion 

6.1. A border line CTTS but with a massive disc 

IM Lupi displays only a modest amount of emission-line 
activity, with a Ha equi valent width which is known to 
vary from 7.5 to 21.5 A (|Batalha fc Basril Il993l ). Studies 
of Ha emission show evidence of accretion, including vari- 
ability: iReipurth et al.1 (|l996h concluded that the Ha fea- 
ture show s an inverse P Cygni p rofile (classification IV- 
Rm) and IWichmann et al.1 (|l999f ) observed a III-R pro- 
file. International Ultraviolet Explorer low dispersion LWP 
spectra show that the Mg II 2 798 A line also y aries, by a 
factor of about 2 in net flux (jValenti et al.|[2003l ) . The rela- 
tively weak emission li nes an d lack of optical vei li ng cau sed 
iFinkenzeller fc Basril l| 19871 ) and iMartin et all (|l994l ) to 
classify this object as a weak-line T Tauri star, although 
our results show it would be more properly categorised as 
a borderline classical T Tauri star. IM Lupi is a relatively 
weak-lined T Tauri star with a large and massive circum- 
stellar disc. Our modelling indicates a large dust disc mass 
of M dust = 10~ 3 M Q , extending up to a radius of 400 AU. 
This large mass is puzzling given the weakness of the Ha 
line, which suggests a low accretion rate. It is however pos- 
sible that the diagnostics of accretion have only been ob- 
served during periods of low or moderate accretion. 

6.2. Disc structure 

Our multi-technique modelling allows us to quantitatively 
constrain most of the geometrical parameters of the disc. 
A flared geometry with a scale height of 10 AU at a ref- 
erence radius of 100 AU is required. The best model has a 
midplane temperature of 14 K at 100 AU. If we assume the 
disc is vertically isothermal (with the temperature equal 
to the midplane temperature), the hydrostatic scale height 
y^kg^T^) /GM^/i (where \i is the mean molecular weight) 
is 12 and 8.5 AU, respectively, for a central mass object of 
0.5 and IMq, respectively. This is in good agreement with 
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Fig. 11. Temperature structure of the disc. The full and 
dashed lines represent the midplane and surface temper- 
atures, respectively. The disc surface at a given radius is 
defined as the altitude above the midplane where the tem- 
perature is maximal. Both temperatures become equal at 
the inner edge of the disc because it is directly heated by 
the stellar radiation. The shaded area represents the tem- 
perature corresponding to the vertical gas scale height of 
the best model, assuming that the disc is vertically isother- 
mal and that the central mass is between 0.5 and 1M Q . 



the scale height deduced from the observations (10 AU), in- 
dicating that the outer parts of the disc are very likely to 
be in hydrostatic equilibrium. Figure ITTI presents the calcu- 
lated midplane temperature compared to the temperature 
corresponding to the gas scale height of the best model. The 
agreement is very good at the inner edge and in the outer 
parts of the disc. In the central parts of the disc (< 30 AU 
and excluding the very inner edge), the midplane tempera- 
ture obtained in our passive disc model is too low to explain 
the scale height required by observations, indicating that an 
additional heating mechanism (like viscous heating) may be 
at play in the disc. 

The inner radius is constrained to be between 0.25 and 
0.40 AU, which corresponds to a maximum dust tempera- 
ture around 1 000 K. The modelling of the SED alone allows 
values of the inner radius down to 0.15 AU corresponding 
to temperatures close to the dust sublimation temperature 
(1500K). The modelling of other observations gives ad- 
ditional constraints on the disc parameters and because of 
the correlations between parameters, this reduces the range 
of possible values for the inner radius. However, scattered 
images and millimetre visibilities constrain the large scale 
structure of the disc. Because we assume that the disc can 
be described in terms of power-laws from the inner edge 
to the outer radius, these constraints affect the derived in- 
ner radius. However, it is very likely that the description of 
the disc using power-laws is too simplistic, and more com- 
plex geometries may slightly shift the inner radius inwards, 
making it compatible with the dust sublimation radius. 

The surface density exponent is found to be close to 
— 1 . This value corresponds to the me dian measurement for 
discs in Taurus- Aurigae obtained by lAndrews &: Williams! 
(2007). This also corresponds to the theoretical value for a 



disc in steady-state accretion ([Hartmann et al.lll998f) . The 
surface density at 5AU is 70g.cm~ 2 . This value is within 
the bro ad peak around the median value of 14g.cm~ 2 for 
Taurus ([Andrews fc Williams! [20071 ). 

Considering a probable stellar mass of sa 1 M Q and a 
gas to dust mass ratio of 100, the disc to star mass ratio is 
ps 0.1, meaning the disc may be unstable through gravita- 
tional collapse. The stability of the disc will depend on the 
surface density. Our derived value of the surface density 
exponent a=-l, as opposed to values of or -2, provides 
disc stabi lity at all radi i according to the Toomre stability 
criterion (lToomrelll964D. Colla pse of gravitationally unsta- 
ble discs (jDurisen et alj [20071 ) is one suggested mode for 
planet formation. The disc of IM Lupi representing about 
1/10 of the star mass, local enhancement of density may 
be sufficient to start planet formation in the disc following 
thi s process. 

Hug hes et al.1 (|2008f) have shown that simultaneous 
studies of the dust continuum and CO emissions in several 
well-studied discs can be reproduced with disc models that 
include a tapered exponential outer edge and not a sharp 
outer radius as we have used here. Current observations 
of IM Lupi do not allow us to study in details the outer 
edge of the disc but the NICMOS image indicates that 
dust grains are still present at radii larger than 400 AU. 
However, as the counter nebulae of the disc is seen in scat- 
tered light images at 0.8 and 1.6/^m, we can get a rough 
estimate of the maximum value for the surface density in 
front of this second nebulae. Dust present at radii larger 
than 400 AU must be optically thin, allowing us to see the 
counter nebulae through the disc and the tentative enve- 
lope. This can be translated into a upper limit for the sur- 
face density of the disc0: S(r « 500 AU) x Ko.S/nm < 1. If we 
assume that the dust composition and grain size distribu- 
tion are the same as in the rest of the disc, this corresponds 
to E(r w 500 AU) < 0.2g.cm" 2 at 500AU, i.e. about 3.5 
times smaller than the extrapolated density from our disc 
model, assuming it extends at radii larger than 400 AU. 
This implies that the density at radii larger than 400 AU 
must decrease significantly faster than the 1/r dependence 
we found for the disc in regions inside 400 AU. 

6.3. Grain growth and dust settling 

As with many other classical T Tauri stars, the slope of IM 
Lupi's millimetre continuum suggests a dust opacity follow- 
ing a law close to re a bs(A) oc A -1 , indicating dust grains in 
the disc are larger than those in the interstellar medium. 
Millimetre photometry mostly traces the thermal emission 
of cold dust in the outer part (> 15 AU) of the disc mid- 
plane (see figure IT2|) , while the hot grains in the central 
regions contribute little to this emission. 

The strong silicate features in the mid-IR spectrum in- 
dicate the presence of micron-sized grains in the disc sur- 
face inside the first central AU. The silicate emission indeed 
comes from the central parts of the disc: 90 % of the emis- 
sion at 10 and 20 /im originates from radii smaller than 1 
and 10 AU, respectively fFigure [T2")) . Larger grains (> 3 /Ltm) 
appear to be almost absent in these regions (section [5. 1| . 

As already mentioned, the slope of the mm continuum 
and the mid-IR silicate features indicate a stratified struc- 



The constraint is stronger at 0.8 than at 1.6 /im, due to 
the larger dust opacity. 
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Fig. 12. Cumulative received fluxes for the best model as 
a function of the distance from the star at a wavelength of 
10 fj,m (full line), 20 /im (dashed line), 70 fj,m (dotted line) 
and 1.3mm (dot-dash line). 



ture for the disc, with large grains in the deeper parts 
of the disc and small grains (< 1 /im) in the surface of 
the disc. This surface layer of small grains remains diffi- 
cult to characterise precisely. We have explored a struc- 
ture where stratification of grain size is vertical and caused 
by settling of larger particles towards the disc midplane. 
Our model with vertical settling allows us to accurately 
reproduce the SED, scattered light images and millime- 
tre emission maps. However, we cannot assess the unique- 
ness of the solution and other explanations may be en- 
visaged. The warm region emitting the 10 micron silicate 
feature is close to the star, while the mm/submm contin- 
uum comes from the whole volume of the disc. We may 
be observing small particles close in and big grains in the 
outer regions of the disc. However, most physical processes: 
grain growth, radial drift by gas drag in the disc mid- 
plane, radiation pressure or stellar wind, will preferentially 
result in larger grains in the central regions and/or re- 
move small grains from these regions. Mid-infr ared inter- 
ferometri c observations (Ivan Boekel et aLll2004l for HAe Be 
stars and iRatzka et al.l 120071 and ISchegerer et al.1 120081 for 
the T Tauri stars TW Hydra and RY Tau respectively) 
have confirmed these trends by showing the presence of 
larger and more processed grains at small spatial scales. 
Timescales for the production of large grains are signif- 
icantly shorter in the central parts of the disc and it is 
difficult to imagine a physical process that will efficiently 
produce large grains in the outer disc without also produc- 
ing them in the inner disc. Even if the current observations 
do not allow us to firmly conclude this point, vertical grain 
size stratification, probably resulting from a combination of 
settling and enhanced grain growth close to the midplane, 
appears to be a more natural explanation, and is, for now, 
our preferred model for the disc of IM Lupi. 

Following the previously described argument, the pres- 
ence of millimetre grains at large radii strongly suggests 
that such large grains are also present in the central AU, 
where the process of grain growth should be more efficient. 
The 10 /im features, dominated by the emission from grains 
of size around 1.5 /im, indicate that the mixing in the disc is 



sufficient to maintain micrometric grains in the disc surface. 
Moreover, the absence of 3 p,m grains in the IRS spectrum 
indicates that the decoupling between gas and dust starts 
for a grain size between 1 and a few /xm, i.e. grains larger 
than this threshold are settled below the surface ri 0Aim = 1. 
Given the high densities in the inner regions of the disc, an 
increase of 1 in optical depth is obtained over a very small 
spatial scale. Even a moderate amount of dust settling is 
sufficient to produce the observed effect. The very low val- 
ues we derive for the settling index confirm that the settling 
remains efficiently counterbalanced by vertical mixing, due 
to turbulence for instance. 

It is very likely that the process of dust settling evolves 
with different timescales as a function of the distance from 
the star. The silicate emission bands provide strong indi- 
cations of the efficiency of dust settling as a mean of re- 
moving grains larger than a few microns from the upper 
layers in the central parts of the disc. The SED also gives 
some insights on the presence of settling in the outer parts. 
The MIPS far-IR fluxes, which are low compared to the 
mid-IR emission, indicate that the disc intercepts a rela- 
tively small fraction of the stellar radiation at large radii 
(> 10 AU) compared to the fraction intercepted at radii 
< 10 AU probed in the mid-IR (see Fig. [T2J). This is ex- 
pected in the case of dust settling (see for instance Fig. 7 
in lD04f ). The decreasing SED of IM Lupi in the mid- infrared 
is ve ry reminiscent of the models including dust settling of 
iDOl Thus, we tried to fit the SED without taking into ac- 
count the mid- infrared fluxes (between 5 and 35 /im), and 
hence not considering the silicate features. The resulting 
Bayesian probabilities still exclude models without dust 
settling, which do not manage to reproduce simultaneously 
both the millimetre and far-infrared fluxes. We conclude 
that dust settling is very likely to occur even in the outer 
reg ions of the disc. 

|D04| predicts that settled discs may become unde- 
tectable in scattered light due to the formation of a self- 
shadowed opacity structure in the outer disc . Thi s is clearly 
not the case for IM Lupi. As noted by |D04| . the self- 
shadowed structure only appears for low values of the tur- 
bulence. This may indicate that the turbulence level in 
IM Lupi is large enough to prevent the disc from self- 
shadowing or that we are observing the outer disc at a 
relatively early stage of the dust settling process. 

6.4. Evolutionary state of IM Lupi & comparison with other 
classical T Tauri stars 

Our modelling has lead to a very detailed picture of the disc 
surrounding IM Lupi. We have already seen that the disc 
structure of IM Lupi is similar to those of other classical 
T Tauri stars. In this section, we compare the signature of 
dust evolution in the disc with the results obtained for other 
T Tauri discs, in order to determine whether it is a singular 
object or whether it can be considered as representative of 
other T Tauri stars. 

Our conclusions on the degree of dust settling in 
IM Lupi result from the strong silicate emission feature 
and shallow millimetre spectral slope. Other objects show 
similar characteristics. Thus, Fig. [13] presents the strength 
of the 10 (im silicate as a function of the exponent of the 
opacity law in the millimetre for a set of T Tauri stars 
(Table [6J . IM Lupi is located in the bulk of T Tauri stars 
and our results can probably be extrapolated to other 
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Fig. 13. Strength of the 10 fim silicate features as a func- 
tion of the millimetre opacity slope for the T Tauri stars 
listed in Table [6l IM Lupi is represented by the red star. 
The effect of grain growth and dust settling are schematised 
by arrows. 

Table 6. Millimetre opacity slopes and strengths of the 
lOjum silicate feature s of T T auri stars. Reference s: (1) = 
ILommen etldl (120071) . (2) = iKessler-Silacci et~aH (120061), 
(3) = iPrzvgodda et all (12003D. (4) = iRodmann et all 
|2006t ). ( 5) = IFurlan et al.l (l2006l) . The source T Cha pre- 
sented in lLommen et all (|2007l ) was not selected, because o f 
the presence of PAH emission. Both lRodmann et all (2006) 
and lLommen et ahl (|2007f) take a contribution of optically 
thick emission into account in their derivation of /3 mm . 



Name 




SpeaktlO/im) 


Ref. 


IM Lup 


0.8 ± 0.25 


1.60 


this work 


HT Lup 


0.4 ± 0.5 


1.39 


(1, 2) 


GW Lup 


0.5 ± 0.5 


1.48 


(1, 2) 


CR Cha 


1.5 ± 0.6 


2.52 


(1, 3) 


WW Cha 


0.8 ± 0.8 


1.92 


(1, 3) 


RU Lup 


0.8 ± 0.5 


1.54 


(1, 2) 


RYTau 


0.8 ± 0.1 


2.75 


(4, 5) 


FT Tau 


0.9 ± 0.3 


1.74 


(4, 5) 


DG Tau 


0.7 ± 0.1 


1.52 


(4, 5) 


UZ Tau E 


0.8 ± 0.1 


1.65 


(4, 5) 


DL Tau 


1.0 ± 0.2 


1.1 


(4, 5) 


CI Tau 


1.3 ± 0.4 


1.5 


(4, 5) 


DO Tau 


0.5 ± 0.1 


1.19 


(4, 5) 


GM Aur 


1.6 ± 0.2 


2.54 


(4, 5) 



sources: detailed analyses of the individual sources are re- 
quired, but it is likely that at least some of them are un- 
dergoing so me dust sett l ing. A semb lable conclusion was 
reache d by iFurlan et all (|2005l l2006f ) and iD'Alessio et~aTl 
(2006| based on similar evidence i.e., the presence in the 
SEDs of 10 and 18 /xm emission silicate bands and the slope 
and fluxes at (sub-)millimetre wavelengths. 

The tentative correlation between the strength of the 
silicate feature and m illimetre spectral index suggested by 
ILommen et al.l (|2007l ) does not appear as clear in our larger 
sample of T Tauri stars. This correlation was interpreted as 
an indication of fast grain growth in both central and outer 
regions of the disc. Indeed, as grains grow from sub-micron 



sizes to several microns, the 10 [im feature becomes weaker 
and less peaked. When they reach millimetre to centimetre 
sizes, the slope of the millimetre emission becomes shal- 
lower. But, as we have shown, the strength of the silicate 
features is strongly related to the degree of dust settling. 
The stronger the settling, the smaller the apparent (i.e. 
probed in the infrared) grain size, making it difficult to ob- 
tain a precise estimate of the actual grain sizes present in 
the midplane. 

The detailed analysis of the silicate features indicates 
a small degree of crystallisation (< 10%) in spite of a 
high mass fr action of micrometric (hence evolved) grains. 
As shown bv lSchegerer et ail (|2006l . see their Fig. 9 for in- 
stance), this is the case for several objects a nd IM Lupi is 
not unique on that aspect. ISchegerer et alJ (|2006h did not 
find any correlation between the degree of crystallisation 
and amount of micron-sized grains. 

IFurlan et al.1 (|2005l . 120061) and ID'Alessio etail (|2006l) 
have estimated the degree of dust settling based on the 
colours and SEDs of a large population of Classical T Tauri 
stars in the Taurus molecular cloud. They claim that to 
make a synthetic SED consistent with the median SED of 
class II objects in Taurus, the dust to gas mass ratio in 
the disc atmosphere should be around 1% (between 0.1 
and 10%) of the ISM ratio. In our modelling, the grain 
size distribution and the dust to gas ratio are continuous 
functions of the height above the midplane. With £ = 0.05, 
the dust to gas ratio starts at 1.5 times the ISM value in 
the disc midplane and reaches 10%, 1% and 0.1% of the 
ISM value at 2, 3 and 6 scale heights, respectively. These 
regions ro ughly corresp o nd to what is defin ed as the disc 
surfac e in IFurlan et all (|2005l . I2006T ) and ID'Alessio et~aT1 
(2006). Although not directly comparable, our estimation 
of the degree of dust settling appears consistent with the 
calculations of these authors. 



7. Summary 

We have constructed a high quality data set of the circum- 
stellar disc of IM Lupi, spanning a wide range of wave- 
lengths, from the optical to the millimetre. All of these 
observations can be interpreted in the framework of sin- 
gle model. A Bayesian analysis of a large grid of models 
allows us to establish strong constraints on all the param- 
eters of the model. Although each individual observation 
provides valuable information on the disc, the presented 
method clearly illustrates the need for multi-wavelength 
and multi-technique studies in order to obtain finer under- 
standing of protoplanetary discs. The conclusions of this 
work are: 

1. The disc structure is very well constrained. The disc ex- 
tends from an inner radius < 0.4 AU, compatible with 
the dust sublimation radius, up to 400 AU. The scale 
height is 10 AU at 100 AU and varies with a flaring index 
of 1.15, indicating that the disc is in hydrostatic equilib- 
rium in its outer parts. The slope of the surface density 
is confined to values close to —1, which is in agr e ement 
with the median value of I Andrews fc Williams! (|2007l) 
and predictions of steady-state accretion disc models. 

2. The millimetre spectral index indicates that grains have 
grown up to a few millimetre in sizes in the disc mid- 
plane. 
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3. The strong silicate emission bands, probing the surface 
of the disc in the central few AUs, also present signa- 
tures of dust evolution. They are dominated by grains 
around 1.5 but are almost devoid of grains larger 
than 3 /xm. 

4. We conclude that a spatial stratification of the dust 
grains, depending on their size is present in the disc. The 
disc of IM Lupi has probably entered the first phases 
of planetary formation: dust grains with sizes several 
orders of magnitude larger than interstellar grains are 
present in the disc and dust settling is probably occur- 
ring, at least in the central parts but potentially also 
in the outer regions of the disc. Models without dust 
settling are excluded but the settling is constrained to 
a low value, suggesting that mixing, by turbulence for 
instance, remains efficient in the disc. 

5. IM Lupi presents signatures of crystalline grains but 
with a low overall degree of crystall inity (< 10%). 
This i s in agreement with the results of ISchegerer et ahl 
(|2006h who found that grain growth and crystallisation 
occurs simultaneously in T Tauri discs, although grain 
growth is dominant. 

6. The simultaneous analysis of the 0.6 and 1.6 //m scat- 
tered light images suggests that light scattering can be 
more forward throwing at short wavelengths. Although 
the data are marginally compatible with compact sili- 
cate grains, this could indicate the presence of grains 
with low refractive indices. They may be fluffy aggre- 
gates (porous grains) and/or the result of the formation 
of ice mantles around grains. Both phenomena are ex- 
pected to occur in discs and additional information (like 
polarisation) is required to distinguish between them. 

The combination of our observations and modelling 
makes IM Lupi one of the best studied protoplanetary disc 
surrounding a solar mass star. With the exception of its low 
accretion signatures, all of the observations indicate that 
IM Lupi is a typical classical T Tauri star. Although signif- 
icant differences are expected between individual objects, 
IM Lupi can probably be considered as good prototype of 
protoplanetary discs for further studies. 

Acknowledgements. Authors would like to thank T. Hill for valuable 
comments on the manuscript. Computations presented in this pa- 
per were performed at the Service Commun de Calcul Intensif de 
l'Observatoire de Grenoble (SCCI) and on the University of Exeter's 
SGI Altix ICE 8200 supercomputer. C. Pinte acknowledges the fund- 
ing from the European Commission's Seventh Framework Program 
as a Marie Curie Intra-European Fellow (PIEF-GA-2008-220891). The 
authors thank the Programme National de Physique Stellaire (PNPS) 
and J'Action SpeciGque en Simulations Numcriques pour V Astronomic 
(ASSNA) of CNRS/INSU, France and Agence Nationale pour la 
Recherche (ANR) of France under contract ANR-07-BLAN-0221, 
for supporting part of this research. This investigation was based, 
in part, on observations made with the NASA/ESA Hubble Space 
Telescope, obtained at the Space Telescope Science Institute (STScI), 
which is operated by the Association of Universities for Research in 
Astronomy, Inc., under NASA contract NAS 5-26555. These observa- 
tions are associated with programs GO/7387 and GO/10177. Support 
for these programs was provided by NASA through grants from STScI. 
This research has made use of the SIMBAD database, operated at 
CDS, Strasbourg, France, and data from the Two Micron All Sky 
Survey (U. Mass, IPAC/CIT) funded by NASA and NSF. Support 
for this work, part of the Spitzer Postdoctoral Fellowship Program, 
was provided by NASA through contracts 1224608, 1230779 and 
1256316, issued by the Jet Propulsion Laboratory, California Institute 
of Technology, under NASA contract 1407. Astrochemistry in Leiden 
is supported by a NWO Spinoza grant and a NOVA grant, and by 
the European Research Training Network "The Origin of Planetary 
Systems" (PLANETS, contract number HPRN-CT-2002-00308). 



References 

Allard, F., Hauschildt, P. H., Alexander, D. R., & Starrfield, S. 1997, 

ARA&A, 35, 137 
Andrews, S. M. & Williams, J. P. 2007, ApJ, 659, 705 
Apai, D., Pascucci, I., Bouwman, J., et al. 2005, Science, 310, 834 
Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 

337, 403 

Barriere-Fouchet, L., Gonzalez, J.-F., Murray, J. R., Humble, R. J., 

& Maddison, S. T. 2005, A&A, 443, 185 
Batalha, C. C. & Basri, G. 1993, ApJ, 412, 363 

Batalha, C. C, Quast, G. R., Torres, C. A. O., et al. 1998, A&AS, 
128, 561 

Beckwith, S. V. W., Henning, T., & Nakagawa, Y. 2000, Protostars 

and Planets IV, 533 
Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, 

AJ, 99, 924 

Bjorkman, J. E. & Wood, K. 2001, ApJ, 554, 615 
Blum, J., Wurm, G., Kempf, S., et al. 2000, Physical Review Letters, 
85, 2426 

Bouwman, J., Henning, T., Hillenbrand, L. A., et al. 2008, ArXiv 
e-prints, 802 

Bouwman, J., Meeus, G., de Koter, A., et al. 2001, A&A, 375, 950 
Bouy, H., Huelamo, N., Pinte, C, et al. 2008, A&A, submitted 
Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 
Bruggcman, D. 1935, Ann. Phys., 24, 636 

Chiang, E. I., Joung, M. K., Creech-Eakman, M. J., et al. 2001, ApJ, 
547, 1077 

D'Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321 
D'Alessio, P., Calvet, N., Hartmann, L., Franco-Hernandez, R., & 

Servm, H. 2006, ApJ, 638, 314 
Dominik, C, Blum, J., Cuzzi, J. N., & Wurm, G. 2007, in Protostars 

and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 783-800 
Dorschner, J., Begemann, B., Henning, T., Jaeger, C, & Mutschke, 

H. 1995, A&A, 300, 503 
Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 
Dullemond, C. P. & Dominik, C. 2004, A&A, 421, 1075 
Dullemond, C. P. & Dominik, C. 2005, A&A, 434, 971 
Durisen, R. H., Boss, A. P., Mayer, L., et al. 2007, in Protostars and 

Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 607-622 
Dutrey, A., Guilloteau, S., & Ho, P. 2007, in Protostars and Planets 

V, ed. B. Reipurth, D. Jewitt, & K. Keil, 495-506 
Evans, N. J., Allen, L. E., Blake, G. A., et al. 2003, PASP, 115, 965 
Evans, N. J., Harvey, M. M., Huard, T. L., et al. 2007, 

Final Delivery of Data from the c2d Legacy Project: IRAC 

and MIPS, Tech. rep., Pasadena: Spitzer Science Center, 

http: / / data.spitzer.caltech.edu/popular / c2d 
Finkenzeller, U. & Basri, G. 1987, ApJ, 318, 823 

Fitzgerald, M. P., Kalas, P. G., Duchene, G., Pinte, C, & Graham, 

J. R. 2007, ApJ, 670, 536 
Fromang, S. & Papaloizou, J. 2006, A&A, 452, 751 
Furlan, E., Calvet, N., D'Alessio, P., et al. 2005, ApJ, 628, L65 
Furlan, E., Hartmann, L., Calvet, N., et al. 2006, ApJS, 165, 568 
Gahm, G. F., Gullbring, E., Fischerstrom, C, Lindroos, K. P., & 

Loden, K. 1993, A&AS, 100, 371 
Ghez, A. M., McCarthy, D. W., Patience, J. L., & Beck, T. L. 1997, 

ApJ, 481, 378 
Gilliland, R. L. 1994, ApJ, 435, L63 
Goldreich, P. & Ward, W. R. 1973, ApJ, 183, 1051 
Graham, J. R., Kalas, P. G., & Matthews, B. C. 2007, ApJ, 654, 595 
Hanner, M. S., Brooke, T. Y., & Tokunaga, A. T. 1998, ApJ, 502, 871 
Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 

495, 385 

Harvey, P., Men'n, B., Huard, T. L., et al. 2007, ApJ, 663, 1149 
Harvey, P. M., Chapman, N., Lai, S.-P., et al. 2006, ApJ, 644, 307 
Henning, T. & Mutschke, H. 1997, A&A, 327, 743 
Hughes, A. M., Wilner, D. J., Qi, C, & Hogerheijde, M. R. 2008, 

ArXiv e-prints, 801 
Hughes, J., Hartigan, P., & Clampitt, L. 1993, AJ, 105, 571 
Hughes, J., Hartigan, P., Krautter, J., & Kelemen, J. 1994, AJ, 108, 

1071 

Irvine, W. M. & Pollack, J. B. 1968, Icarus, 8, 324 
Jaeger, C, Molster, F. J., Dorschner, J., et al. 1998, A&A, 339, 904 
Johansen, A., Henning, T., & Klahr, H. 2006, ApJ, 643, 1219 
Johansen, A., Oishi, J. S., Low, M.-M. M., et al. 2007, Nature, 448, 
1022 

Jones, A. P., Tielens, A. G. G. M., & Hollenbach, D. J. 1996, ApJ, 
469, 740 



20 



C. Pinte et al.: Probing dust grain evolution in IM Lupi's circumstellar disc 



Kemper, F., Vriend, W. J., & Tielens, A. G. G. M. 2004, ApJ, 609, 
826 

Kessler-Silacci, J., Augereau, J.-C., Dullcmond, C. P., ct al. 2006, 
ApJ, 639, 275 

Koike, C, Tsuchiyama, A., Shibai, H., et al. 2000, A&A, 363, 1115 
Krautter, J. 1992, The Star Forming Region in Lupus (Low Mass Star 

Formation in Southern Molecular Clouds), 127 
Krist, J. E., Burrows, C. J., Stapelfeldt, K. R., & WFPC2 Id Team. 

1997, Bulletin of the American Astronomical Society, 29, 1215 
Krist, J. E. & Hook, R. N. 1997, in The 1997 HST Calibration 

Workshop with a New Generation of Instruments, p. 192, ed. 

S. Casertano, R. Jedrzcjewski, T. Keyes, & M. Stevens, 192 
Lahuis, F., Kessler-Silacci, J. E., Evans, N. J., I., et al. 2006, 

c2d Spectroscopy Explanatory Supplement, Tech. rep., Pasadena: 

Spitzer Science Center 
Lay, O. P., Carlstrom, J. E., & Hills, R. E. 1997, ApJ, 489, 917 
Lissauer, J. J. & Stevenson, D. J. 2007, in Protostars and Planets V, 

cd. B. Reipurth, D. Jewitt, & K. Keil, 591-606 
Lommcn, D., Wright, C. M., Maddison, S. T., et al. 2007, A&A, 462, 

211 

Lucy, L. B. 1999, A&A, 345, 211 

Martin, E. L., Rebolo, R., Magazzu, A., & Pavlenko, Y. V. 1994, 

A&A, 282, 503 
Mathis, J. S. & Whiffen, G. 1989, ApJ, 341, 808 

Meri'n, B., Augereau, J.-C, van Dishoeck, E. F., et al. 2007, ApJ, 661, 
361 

Min, M., Dominik, C, Hovenier, J. W., de Koter, A., & Waters, 

L. B. F. M. 2006, A&A, 445, 1005 
Molster, F. J., Waters, L. B. F. M., & Tielens, A. G. G. M. 2002, 

A&A, 382, 222 

Natta, A., Testi, L., Calvet, N., et al. 2007, in Protostars and Planets 

V, ed. B. Reipurth, D. Jewitt, & K. Keil, 767-781 
Nuernberger, D., Chini, R., & Zinnecker, H. 1997, A&A, 324, 1036 
Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 
215 

Padgett, D. L., Cicza, L., Stapelfeldt, K. R., et al. 2006, accepted for 

publication in ApJ, astro-ph/0603370 
Pinte, C, Fouchet, L., Menard, F., Gonzalez, J.-F., & Duchcne, G. 

2007, A&A, 469, 963 
Pinte, C, Menard, F., Duchene, G., & Bastien, P. 2006, A&A, 459, 

797 

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, 
B. P. 1992, Numerical recipes. The art of scientific computing 
(Cambridge: University Press, 1992, 2nd cd.) 

Przygodda, F., van Boekel, R., Abraham, P., et al. 2003, A&A, 412, 
L43 

Ratzka, T., Leinert, C, Henning, T., et al. 2007, A&A, 471, 173 
Reipurth, B., Pedrosa, A., & Lago, M. T. V. T. 1996, A&AS, 120, 229 
Rodmann, J., Henning, T., Chandler, C. J., Mundy, L. G., & Wilncr, 

D. J. 2006, A&A, 446, 211 
Safronov, V. S. & Zvjagina, E. V. 1969, Icarus, 10, 109 
Schcgcrer, A., Wolf, S., Voshchinnikov, N. V., Przygodda, F., & 

Kessler-Silacci, J. E. 2006, A&A, 456, 535 
Schegerer, A. A., Wolf, S., Ratzka, T., & Leinert, C. 2008, A&A, 478, 

779 

Schneider, G., Silverstone, M. D., & Hines, D. C. 2005, ApJ, 629, L117 
Schneider, G., Silverstone, M. D., Hines, D. C, et al. 2006, ApJ, 650, 
414 

Schneider, G., Wood, K., Silverstone, M. D., et al. 2003, AJ, 125, 1467 
Schraplcr, R. & Henning, T. 2004, ApJ, 614, 960 
Servoin, J. & Piriou, B. 1973, Phys. Stat. Sol. (b), 55, 677 
Tachihara, K., Dobashi, K., Mizuno, A., Ogawa, H., & Fukui, Y. 1996, 
PASJ, 48, 489 

Thamm, E., Stcinacker, J., & Henning, T. 1994, A&A, 287, 493 
Toomre, A. 1964, ApJ, 139, 1217 

Valenti, J. A., Fallon, A. A., & Johns-Krull, C. M. 2003, ApJS, 147, 
305 

van Boekel, R., Min, M., Leinert, C, et al. 2004, Nature, 432, 479 
van Boekel, R., Min, M., Waters, L. B. F. M., et al. 2005, A&A, 437, 
189 

van Kcmpen, T. A., van Dishoeck, E. F., Brinch, C, & Hogerheijde, 

M. R. 2007, A&A, 461, 983 
Voshchinnikov, N. V. & Henning, T. 2008, A&A, 483, L9 
Voshchinnikov, N. V., Il'in, V. B., & Henning, T. 2005, A&A, 429, 

371 

Watson, A. M., Stapelfeldt, K. R., Wood, K., & Menard, F. 2007, in 
Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 
523-538 



Weidenschilling, S. J. 1977, MNRAS, 180, 57 

Wichmann, R., Bastian, U., Krautter, J., Jankovics, I., & Rucinski, 

S. M. 1998, MNRAS, 301, L39 
Wichmann, R., Covino, E., Alcala, J. M., et al. 1999, MNRAS, 307, 

909 

Wright, E. L. 1987, ApJ, 320, 818 



